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1  Introduction 

This  report  contains  two  related  manuscripts  on  the  subject  of  mechanism 
geometry.  The  first  paper  addresses  general  theory,  the  second  addresses 
practical  computations.  Together  they  suggest  a  general  approach  for  ef¬ 
ficient  parametrization  of  spatial  mechanism  configuration  spaces.  Before 
describing  the  contents  of  these  papers  we  discuss  some  of  the  steps  that  led 
us  to  this  point. 

About  ten  years  ago  we  started  developing  a  general  control  architecture 
for  launch  vehicles  as  part  of  the  Adaptive  Guidance,  Navigation  and  Control 
contract.  The  objective  was  to  reduce  the  time  and  cost  of  producing  flight 
software  for  launch  vehicles. 

One  of  the  diflBcult  features  of  laimch  vehicles,  from  a  control  perspective, 
is  their  limited  margins  in  both  performance  and  robustness.  Because  of  lim¬ 
ited  performance  margin,  the  ascent  trajectory  must  be  carefully  tailored  for 
each  flight,  depending  on  payload  and  mission.  Because  of  limited  robustness 
margin,  the  control  laws  must  be  reanalyzed  and  perhaps  changed  depending 
on  payload  and  mission.  In  order  to  reduce  the  time  and  effort  required  for 
the  flight  software  development,  we  decided  to  automate  the  guidance  and 
control  design  tasks  in  a  software-development  environment  that  included  a 
multibody  simulation  of  the  launch  vehicle.  The  multibody  model  reflected 
propellant  slosh,  engine  gimbaling  and  vehicle  flex  dynamics.  For  this  reason 
we  wanted  to  find  an  efiicient  simulation  of  our  multibody  system. 

The  multibody  model  we  adopted  for  the  laimch  vehicle  had  one  feature 
that  made  the  kinematic  problem  especially  easy  -  there  were  no  closed  loops. 
In  such  a  case  the  configuration  space  is  simply  the  product  of  the  individual 
configuration  spaces  of  the  joints  linking  the  bodies  together.  We  were  using 
revolute  and  spherical  joints  only,  so  the  configuration  space  turned  out  to 
be  a  product  of  (closed  subsets  of)  circles  and  orthogonal  groups.  Easy  to 
parametrize.  The  key  technical  insight  for  this  problem  was  the  realization 
that  the  dynamic  computations  could  be  reduced  from  order  to  order  N , 
where  N  is  the  number  of  degrees  of  freedom.  Having  figured  out  how  to  do 
that  for  our  launch  vehicle  simulation,  we  survived  our  first  serious  encounter 
with  multibody  systems. 

The  next  multibody  problem  we  encountered  was  a  true  mechanism  prob¬ 
lem  -  a  six-degree-of-freedom  hand  controller.  The  hand  controller  mecha¬ 
nism  was  a  general  Stuart  platform  with  three  legs  connecting  the  movable 
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top  plate  to  the  base.  The  details  of  the  configuration  are  not  important  for 
our  discussion  here,  what  is  critical  is  that  the  motions  of  the  various  bodies 
maVing  up  the  mechanism  were  constrained  by  loop-closure  conditions.  In 
this  situation  is  it  is  no  longer  true  that  the  configuration  space  is  the  prod¬ 
uct  of  the  individual  configuration  spaces  of  the  joints  coimecting  the  bodies. 
In  fact,  for  general  mechanisms  of  this  type,  the  configuration  spaces  may  be 
very  complicated  and  need  not  be  easy  to  parametrize.  This  state  of  affairs 
can  make  it  difficult  to  design  a  mechanism  for  an  engineering  application, 
because  it  is  not  easy  to  write  down  and  anal3rze  dynamic  equations  without 
first  parametrizing  the  configuration  space  (methods  do  exist,  but  they  are 
not  so  easy  to  use  for  system  design).  In  any  case,  after  some  effort  we  were 
able  to  parametrize  our  hand-controller  mechanism  and  so  we  survived  our 
second  encounter  with  multibody  systems. 

But  our  curiosity  was  piqued.  The  analysis  of  the  hand  controller  took 
longer  than  we  thought  it  should,  and  the  approach  that  finally  worked  was 
imsystematic.  Surely  there  was  a  better  approach,  one  that  we  could  learn  be¬ 
fore  our  third  encounter.  On  looking  into  the  matter,  however,  we  fotmd  that 
no  general,  systematic  approach  was  known.  In  fact,  there  were  single-loop 
one-DOF  mechanisms  (the  7R)  for  which  the  most  authoritative  reference  we 
could  find  (Duffy’s  book)  had  no  solution.  In  the  1970’s,  the  7R  problem  was 
called  the  Mount  Everest  of  mechanism  problems  by  one  of  the  experts  in 
the  field.  As  luck  would  have  it,  the  Chinese  researchers  Lee  and  Liang  had 
just  solved  the  7R  problem  in  1986  (it  was  now  the  early  1990’s).  Ignorant 
of  their  results,  we  set  out  to  develop  a  more  systematic  approach  to  general 
mechanism  problems,  and  we  decided  to  use  the  7R  as  our  acid  test  for  the 
general  method. 

After  a  little  work  we  found  an  approach  that  was  at  the  same  time  sys¬ 
tematic  and  computationally  tractable.  We  worked  out  the  details  for  the 
special  case  of  the  7R  mechanism  and  submitted  our  results  for  publication  in 
June  of  1993.  It  was  then  we  learned  of  the  Chinese  solution  of  the  7R  prob¬ 
lem,  so  we  could  not  claim  primacy  for  the  7R  solution.  But  our  approach 
to  the  problem  was  original  enough  to  warrant  publication,  so  we  revised 
the  manuscript  and  resubmitted  it.  It  appeared  as  “  A  New  Computational 
Algorithm  for  7R  Spatial  Mechanisms  ”  in  Mech.  Mach.  Theory  Vol.  31, 
No.  1,  pp.  23-43,  1996. 

In  the  five  years  since  we  first  applied  our  method  to  the  7R  problem  we 
have  continued  developing  both  the  theory  and  the  computational  methods 
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associated  with  the  general  approach.  The  fundamental  theory  is  described 
in  full  generality  in  the  first  paper  included  in  this  report.  This  first  paper 
is  complete  as  it  stands  -  ready  for  publication.  The  second  paper  is  not 
ready  for  publication,  it  is  the  current  draft  of  a  working  document  on  the 
computational  technique.  All  of  the  methods  discussed  in  the  second  paper 
have  been  tested  on  numerical  examples  with  experimental  software,  so  we 
are  confident  of  useful  results  even  though  we  do  not  yet  know  the  full  range 
of  applicability.  We  expect  to  continue  research  into  this  area  in  the  years  to 
come. 


1.1  Discussion  of  Key  Results 

The  first  paper  is  titled  “  Multi-Affine  Modular  Systems  for  Mechanism  Fam¬ 
ilies.”  A  summary  of  contents  can  be  foimd  in  the  abstract  and  introduction 
of  that  paper.  Here  we  outline  the  key  steps. 

The  first  step  is  to  define  precisely  what  a  mechanism  of  (R,P,S)-t3q)e 
is.  A  minimal  set  of  parameters  needed  to  describe  a  mechanism  of  this 
type  is  then  defined.  We  then  describe  what  is  meant  by  the  configuration 
space  of  a  mechanism,  and  develop  ideas  for  computing  the  dimension  of  the 
configuration  space  firom  the  parameters  describing  the  mechanism. 

Some  illustrative  examples  are  then  developed  to  show  the  difficulty  with 
the  usual  formula  for  the  dimension  of  the  configuration  space  from  crude 
features  of  the  mechanism  data  alone.  This  difficulty  motivates  consideration 
of  families  of  mechanisms,  and  leads  to  the  distinction  between  complete 
and  overconstrained  families.  A  general  theory  of  complete  families  is  then 
developed,  and  it  is  shown  that  the  usual  formula  for  the  dimension  of  the 
configuration  space  applies  to  a  generic  mechanism  in  such  a  family.  It  is 
worth  pointing  out  that  nongeneric  mechanisms  in  complete  families  and  even 
overconstrained  mechanisms  are  often  (approximately)  encountered  in  real- 
world  systems  (e.g.  door-locks,  axles),  so  the  mathematical  ideas  developed 
here  are  of  practical  significance,  even  though  they  are  not  usually  considered. 

Next,  for  a  generic  mechanism  in  a  complete  family  we  describe  the  struc¬ 
ture  of  a  branched  cover  for  the  configuration  space  over  a  base  space  com¬ 
posed  of  a  product  of  joint  configuration  spaces  (this  structure  often  arises, 
but  not  always  -  consider  three  bodies  connected  by  spherical  joints  for  a 
counterexample).  The  physical  significance  of  the  branch  locus  and  the  cov¬ 
ering  degree  of  this  branched  cover  (when  it  exists)  are  discussed.  It  is  an 
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interesting,  open  problem  to  find  a  formula  for  the  covering  degree  directly 
from  the  mechanism  parameters  and  the  joint-data  of  the  base  space.  In 
this  situation,  to  compute  coordinates  on  the  configuration  space,  we  can 
compute  local  sections  of  the  branched  cover  over  the  base  space. 

Whether  or  not  the  branched  cover  over  the  product  of  joints  exists,  we 
find  a  covering  set  of  local  parametrizations.  Using  a  special  set  of 
joint  coordinates,  we  show  how  a  modular  system  of  multiafiine  equations 
describing  the  configuration  space  can  be  generated.  The  generation  of  the 
multiaflfine  equations  is  the  central  theme  of  the  third  section  of  this  paper. 
Methods  for  solving  equations  of  this  type  are  the  topic  of  the  second  paper 
in  this  report. 

The  second  paper  is  titled  “A  Polynomial-Runtime  Algorithm  for  Multi- 
Afiine  Equations.”  The  central  theme  of  this  second  paper  is  a  polynomial¬ 
runtime  ^gorithm  for  solving  the  system  of  equations  developed  in  the  first 
paper.  In  this  case,  the  polynomial  time  condition  means  polynomial  in  the 
number  of  solutions  to  the  system  of  equations.  In  general,  the  number  of 
solutions  of  the  system  of  equations  grows  combinatorially  with  the  number  of 
equations,  so  the  runtime  of  the  algorithm  grows  faster  than  any  polynomial 
function  of  the  number  of  equations.  Given  the  growth  in  the  number  of 
solutions,  however,  there  is  no  avoiding  this  problem  (just  writing  down  the 
solutions  to  the  equations  has  superpolynomial  growth).  So,  an  algorithm 
polynomial  in  the  number  of  solutions  is  the  best  one  can  hope  for.  It  is 
worth  noting  that  Groebner  basis  methods  are  worse  than  polynomial  in  the 
number  of  solutions,  so  the  algorithm  described  here  is  an  improvement  over 
Groebner  basis. 

The  text  is  divided  into  two  chapters.  The  first  chapter  concerns  solu¬ 
tion  methods  when  the  number  of  equations  is  the  same  as  the  number  of 
unknowns.  In  this  case  the  number  of  solutions  for  a  generic  system  is  given 
by  a  known  formula,  and  the  solution  method  uses  that  formula. 

The  two  basic  techniques  are: 

1.  Dilation  -  generating  many  equations  of  higher  degree  from  the  original 
multiafiine  system 

2.  Defiation  -  removing  extraneous  roots  from  the  higher-order  system 

When  these  two  methods  are  applied  to  systems  of  multiafiine  equations, 
the  result  is  a  generalized  eigenproblem  of  size  equal  to  the  number  of  solu- 


tions.  The  eigenproblem  can  be  solved  in  polynomial  time,  hence  the  poly¬ 
nomial  runtime  of  the  overall  algorithm. 

The  second  chapter  concerns  solution  methods  for  multiaffine  systems 
when  the  number  of  equations  exceeds  the  number  of  unknowns.  Such  sys¬ 
tems  are  called  overdetermined.  Systems  of  this  type  are  generated  by  the 
method  described  in  the  first  paper,  hence  our  interest  in  them.  When  coef¬ 
ficients  are  generic,  overdetermined  systems  have  no  solutions,  so  the  coefii- 
cients  that  arise  in  the  mechanism  problem  must  satisfy  some  special  set  of 
relations  that  allows  solutions  to  exist. 

Overdetermined  systems  in  n  variables  tend  to  have  fewer  solutions  than 
systems  with  n  equations  in  n  variables,  so  in  theory  the  solution  algorithm 
could  be  quicker.  In  practice,  the  additional  structure  usually  makes  the  so¬ 
lution  more  diflScult  -  especially  in  those  cases  where  the  number  of  solutions 
is  not  know  a  priori. 

For  the  overdetermined  systems  arising  in  mechanism  .anals^is  we  have 
developed  a  third  technique  described  in  section  2.3.  It  is  an  iterative  method 
for  reducing  the  size  of  a  nonsquare  generalized  eigenproblem  to  a  smaller 
square  one.  A  tentative  name  for  this  procedure  is  ablation,  but  a  better 
name  might  be  found  when  the  theory  is  better  developed. 

The  theory  in  the  case  of  overdetermined  systems  clearly  needs  more 
work.  By  implementing  these  algorithms  in  experimental  code  and  applying 
them  to  examples  we  know  the  methods  work  at  least  for  some  problems, 
but  we  have  yet  to  determine  the  complete  range  of  applicability. 
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Abstract 


This  papa*  describes  part  of  our  six-step,  general  approach  for  mechanism  kinematic  analysis. 
The  details  for  a  particular  application  for  the  7R  spatial  mechanism  are  presented  in  [ME]. 
Here  we  concentrate  on  the  derivation  of  the  modular  system  for  a  general  class  of  spatial 
mechanisms. 

The  quaternion-pair  approach  described  below  can  be  applied  to  spatial  mechanisms  endowed 
with  spherical,  prismatic  and  cylindrical  joints  as  well  as  the  revolute  joints  illustrated  in 
[ME].  Systems  of  multi-affine  equations  are  obtained  in  the  same  way  for  single  and  multi¬ 
loop  mechanisms.  In  these  equations,  the  number  of  independent  variables  is  equal  to  the  sum 
over  all  joints  of  the  joint  degrees  of  freedom,  so  such  a  system  is  minimal  in  that  sense. 
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1.  INTRODUCTION 


The  algebraic  analysis  of  mechanism  kinematics  is  an  old  subject  that  remains  an  active  area 
of  research  today.  Some  of  the  current  theoretical  interest  has  been  motivated  by  the  growing 
power  in  today’s  computers,  to  analyze  examples  that  were  simply  too  complicated  to  exam¬ 
ine  by  hand.  A  case  in  point  is  the  general  1-DOF,  spatial  7R  mechamsm  that  was  recently 
solved  by  Lee  and  Liang  iLeel,2]  in  the  late  1980s.  The  subject  of  this  paper  is  a  new 
approach  to  kinematic  analysis  of  mechanisms  that  addresses  a  special  but  important  class  of 
examples,  including  the  7R  mechanism.  Here  we  describe  our  general  approach  from  a 
mathematical  point  of  view.  The  main  result  is  a  computational  technique  for  generating 
modular  systems  of  equations  [Macaulay]  of  a  very  special  type. 

Before  embarking  on  the  general  theory  we  present  a  simple  example,  to  demonstrate  the 
nature  of  the  the  problem  being  solved. 

Consider  a  4-bar,  planar  linkage.  One  of  the  bars  is  assumed  fixed  (ground),  while  the  other 
three  are  joined  in  pairs  by  revolute  joints  to  form  a  closed  circuit.  For  general  values  of  the 
bar  lengths,  this  simple  mechanism  has  one  degree  of  freedom.  Any  one  of  the  joint  angles 
can  be  used  as  a  local  coordinate  for  the  one-dimensional  configuration  space,  pick  one  of  the 
angles  on  the  base  and  call  it  63.  The  problem  is  to  find  the  values  of  Gq,  0i,  and  02  as  a 
function  of  03. 

There  is  a  standard  technique  for  converting  this  problem  into  a  system  of  algebraic  equations. 
For  each  fixed  value  of  0^  there  is  a  simple  geometric  construction  to  find  the  allowed 
configurations.  The  problem  reduces  to  finding  the  points  of  intersection  of  two  circles  in  the 
plane,  and  fiom  the  geometry  it  is  clear  that  there  are  only  two  configurations  possible  for 
each  value  of  0^.  The  obvious  way  to  compute  the  coordinates  of  these  two  points  (and 
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thereby  determine  the  other  joint  angles)  is  to  write  down  the  pair  of  quadratic  equations  in 
the  coordinates  (x,y)  of  the  intersection  point  and  find  all  solutions.  An  interesting  feature 
arises  from  the  algebraic  viewpoint  —  by  Bezout’s  theorem  we  might  expect  to  find  foiu:  solu¬ 
tions  rathCT  than  two.  This  apparent  inconsistency  is  resolved  by  a  well-known  explanation. 
Two  of  the  four  algebraic  solutions  are  finite,  possibly  corresponding  to  real  geometric  solu¬ 
tions  (if  real  solutions  exist),  while  the  remaining  two  solutions  are  always  meaningless  for 
the  kinematic  problem  ~  they  are  the  circular  points  at  infinity.  This  problem  of  extraneous 
roots  is  easily  handled  in  the  simple  examples,  but  the  analogous  problem  for  more  complex 
tnp^hanismg  has  not  been  solved  in  general.  The  work  by  Duffy  [Duffy  1]  is  the  most 
comprehensive  work  in  which  this  type  of  problem  has  been  addressed,  but  all  the  spatial 
mechanisms  analyzed  there  are  1-DOF,  single-loop  with  prismatic  and  revolute  joints.  The 
difficulties  presented  by  the  7R  problem,  for  example,  can  be  traced  in  the  evolutionary  path 
[Freudenstein],  [Roth],  [Duffy2],  [Lee  1,2]  that  led  to  the  correct  solution  in  this  case. 

The  solution  presented  in  [Lee  1,2]  for  the  7R  problem  was  not  very  systematic,  leaving  open 
the  question  of  how  to  proceed  for  more  complicated  (e.g.  multiloop)  spatial  algorithm.  To 
address  the  problem  of  general  multiloop  spatial  mechanisms,  we  have  developed  the  alterna¬ 
tive  approach  presented  in  this  paper.  This  method  was  first  applied  to  the  7R  mechanism  in 
[ME].  The  general  mathematical  theory  is  presented  here. 


To  give  a  flavor  of  our  new  approach,  let  us  return  to  the  four-bar  example.  By  the  methods 
described  in  Section  3,  we  derive  systems  of  multi-affine  equations  in  the  four  variables 


Zj,  Z2,  Z3,  Z4  where 


Zj  =  tan(-^) 


(1.1) 


In  the  special  case  of  planar  4-bars,  we  obtain  three  equations  for  each  of  the  four  angle 
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parameters  —  obviously  these  are  not  independent.  It  turns  out  that  only  two  of  these  three 
(the  translation  equations)  represent  real  constraints  for  the  parameters  Zj.  To  show  the  sym¬ 
metry  of  the  construction  process  we  use  a  variable  index  notation  z^.^2>  ^k+i» 
sider  the  subscripts  in  the  integers  mod  4  for  each  integer  value  of  k. 


For  each  fixed  value  of  7^+3,  we  can  solve  for  Zjj+2.  Zk+i’  ^k-  Let  the  leg  lengths  be  Ij,  I2, 
I3,  and  I4.  Two  components  of  the  translation  equations  (starting  at  joint  k)  are: 

=  0eR2  (1.2) 


A(Zij^3)  +  Zi^+2  B(Z]j+3) 


1 

^k+1 


which  is  a  generalized  eigenvalue  problem  with: 


A(Zk+3)  = 


W+lk+3+W+lk+l 

Zk+3(”lk+4+lk+3+W2+lk+l) 


2^k+3Hk+4+W3'*'lk+2~Wl) 

“W“lk+3~lk+2+lk+l 


(1.3) 


B(Zk+3)  = 


Zk+3(“lk+4+lk+3“lk+2~lk+l) 

~lk+4“lk+3+lk+2+lk+l 


~lk44“lk+3+lk+2~lk+l 

Zk+3(lk44“lk+3+lk+2“lk+l) 


(1.4) 


To  obtain  the  remaining  variable,  zjj,  we  can  use  an  additional  affine  equation  obtained  by 
starting  the  translation  at  another  joint.  By  starting  at  each  of  the  4  joints,  the  translation 
equations  give  a  set  of  two  scalar  multi-affine  equations  for  each  choice  of  k=l, 2,3,4  for  a 
total  of  8  equations.  Since  the  k*  equations  do  not  involve  z^,  we  can  obtain  8  more  multi- 
affine  equations  by  multiplying  by  zjj,  for  a  total  of  16  multiaffine  equations.  This  set  of 
multi-affine  equations  is  an  example  of  our  main  result. 


In  any  application  there  is  more  work  to  be  done  once  the  system  of  equations  has  been  con¬ 
structed.  Let  us  consider  the  last  step  for  the  planar  four-bar.  From  the  16  multi-affine 
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equations  there  are  many  ways  to  select  a  set  of  four  independent  ones  to  obtain  a  4  x  4  gen¬ 
eralized  eigenvalue  problem  involving  all  the  Zj. 


(1.5) 


The  determinant  of  the  4  x  4  matrix  on  the  left-hand  side  of  equation  (1.5)  must  vanish,  giv¬ 
ing  us  a  single  degree-4  polynomial  relation  between  Zj^  and  There  are  several  ways  to 
proceed  fix)m  here. 

First,  by  proper  choice  of  multi-affine  equations  in  constructing  the  matrix,  one  can  find  spe¬ 
cial  structure  in  the  matrices  L  and  N  that  leads  to  degeneracies  in  the  determinantal  polyno¬ 
mial.  The  determinant  becomes  the  square  of  a  degree-two  polynomial  so  only  two  distinct 
solutions  exist. 

Alternately,  one  can  take  the  point  of  view  that  there  are  many  such  determinantal  polynomi¬ 
als  in  the  samp,  two  variables.  The  polynomial  we  want  is  the  greatest  common  divisor  of  the 
complete  set. 

In  any  case,  when  all  the  multiaffine  equations  are  considered  together,  it  is  found  that  only 
two  solutions  exist 

The  reader  is  invited  to  try  his  hand  in  exploring  this  simple  case  computationally.  The  gen¬ 
eral  situation  for  more  complicated  mechanisms  remains  an  interesting  research  problem. 
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2.  PROBLEM  STATEMENT  AND  NOTATION 


We  will  be  considering  spatial  mechanisms  consisting  of  rigid  bodies  connected  by  joints  of 
the  following  type: 

1)  revolute  (type  R) 

2)  prismatic  (type  P) 

3)  spherical  (type  S) 


We  assume  the  reader  is  familiar  with  the  above  terminology  -  the  mathematical  characteriza¬ 
tion  of  these  joints  is  discussed  in  the  next  paragraph.  The  standard  enginnering  perspective  is 
presented  in  [Duffy].  It  should  be  observed  that  some  other  joint  types  (e.g.  cylindric),  some¬ 
times  considered  separately  by  engineers,  can  be  constructed  from  the  basic  joint  types  above. 

The  common  feature  shared  by  the  above  three  joint  types  can  be  stated  succinctly:  the  rela¬ 
tive  motions  allowed  by  these  joints  form  real-algebraic  subgroups  of  the  group  E(3)  of 
Euclidean  motions  in  R^.  The  relative  motions  allowed  by  an  R-joint  are  those  that  fix  a  point 
in  space  and  a  line  containing  that  point.  The  relative  motions  allowed  by  a  P-joint  form  a  1- 
parameter  translation  group.  The  relative  motions  allowed  by  an  S-joint  are  those  that  fix  a 
point  in  space. 

2.1  Mechanism  Definition 

To  specify  a  mechanism  completely  (for  our  purposes),  we  need  the  following  data: 

1)  A  set  W  of  rigid  bodies 

2)  A  subset  J  of  W  x  W;  points  in  J  are  pairs  of  bodies  that  are  connected 
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3)  A  joint  (R,  P,  or  S)  for  each  point  in  J 

4)  Data  specifying  the  relative  constraint  imposed  by  each  joint  on  its  body-pair 
The  number  of  bodies  (cardinality  of  W)  is  denoted  N. 

We  e3q)lain  "mechanism  specification"  by  a  practical  analogy.  Imagine  a  box,  full  of  parts,  on 
which  is  written  "some  assembly  retjuiied."  Included  in  the  box  is  an  instruction  booklet  that 
tells  how  to  identify  each  part  umquely  (identify  the  points  in  W).  The  booklet  also  tells 
which  part  should  be  attached  to  another  part  (identify  the  points  in  J),  and  the  means  of 
attachment  is  specified  as  well  (joint  type  for  each  point  in  J).  On  each  body,  the  manufac¬ 
turer  has  thoughtfully  provided  fixtures  appropriate  for  the  specified  joints.  The  data  specify¬ 
ing  the  relative  constraint  imposed  by  each  joint  is  exactly  the  information  required  by  the 
machine  operator  who  endowed  each  body  with  its  joint  fixtures.  The  reader  can  imagine 
assembling  the  mechanism,  one  joint  at  a  time  (perhaps  according  to  some  recommended 
order),  until  all  the  bodies  are  connected  by  the  appropriate  joints. 

Anyone  who  has  bought  a  box  full  of  parts  on  which  is  written  "some  assembly  required"  will 
no  doubt  be  wondering  if  the  parts  in  our  hypothetical  box  really  do  fit  together.  In  answer  to 
that  question  we,  like  most  manufacturers,  merely  give  assurance  that  the  parts  were  designed 
to  fit  together.  We  leave  some  constraction  details  to  the  customer  who  should  find  that,  for 
some  configurations  of  pieces  during  assembly,  the  parts  really  do  fit  together  as  advertised 
(provided  the  dimensions  of  the  pieces  are  within  their  design  tolerances). 

2.2  The  Configuration  Space  of  a  Mechanism 

Once  the  mechanism  is  assembled,  the  individual  bodies  within  it  may  be  moved  continuously 
through  a  range  of  distinct  positions  and  orientations  relative  to  one  another.  When  we  con¬ 
sider  motions,  one  specified  body  is  assumed  fixed  in  space  (attached  to  ground)  to  eliminate 
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motions  of  the  entire  mechanism  as  a  rigid  body.  The  set  C  of  all  possible  distinct,  relative 
positions  of  the  bodies  in  the  assembled  mechanism  is  called  the  configuration  space.  Let  us 
clarify  what  we  mean  by  the  set  of  possible  distinct,  relative  positions  of  the  bodies. 

We  return  to  our  practical  analogy.  Unlike  most  manufacturers,  we  explicitly  state  that  the 
pieces  might  not  fit  together  if  their  dimensions  are  outside  design  tolerances.  This  can  hap¬ 
pen,  for  example,  in  a  4-bar  planar  linkage  where  the  length  of  one  leg  is  greater  than  the  sum 
of  the  lengths  of  the  other  three.  A  more  subtle  concern  is  that  different  (i.e.  nonisotopic) 
embeddings  of  the  mechanism  in  could  result  if  two  identical  kits  are  assembled  by  pro¬ 
cedures  involving  different  intermediate  configurations.  Moving  the  mechanism  too  far  in  one 
direction  will  cause  the  bodies  to  collide.  Problems  of  this  sort  can  be  very  complicated  if 
considered  in  the  context  of  global  embeddings  --  the  topological  questions  alone  lead  to 
problems  in  knot  and  braid  theory.  We  look  instead  at  the  simpler  problem  in  which  intersec¬ 
tions  are  allowed  to  occur  and  the  various  bodies  are  allowed  to  pass  through  each  other. 
That  is,  we  allow  mechanism  configurations  into  the  set  C  even  if  they  are  immersed  but  not 
necessarily  embedded.  In  addition,  we  allow  some  complex  solutions.  Allowing  complex  solu¬ 
tions  (at  least  for  some  values  of  the  part  dimensions)  is  an  artifice  mandated  by  our  computa¬ 
tional  techniques.  The  set  of  real,  embedded  (i.e.  physically  realizable)  configurations  is  a 
subset  of  the  configuration  space  C  we  have  defined  here. 

2.3  Families  of  Mechanisms 


We  now  define  the  concept  of  a  family  of  mechanisms. 

Definition:  By  a  family  of  mechaitisms,  we  mean  a  collection  of  mechanisms  having  the  same 
number  of  bodies,  isomorphic  sets  J,  and  the  same  joint  type  for  the  corresponding  points  in  J. 

In  other  words,  the  only  differences  between  two  mechanisms  in  the  same  family  are: 
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1)  sizes  and  shapes  of  the  N  bodies  (independent  of  joint  data) 

2)  data  involving  fixed-point  location  and  line  orientation  at  each  R-joint 

3)  direction  of  translation  for  each  P-joint 

4)  fixed-point  location  for  each  S-joint 

Because  we  are  only  concerned  with  immersed  mechanisms,  we  really  do  not  need  to  consider 
the  exact  sizes  and  shapes  of  the  bodies  (parameters  of  the  first  type).  Each  body  might  as 
well  be  all  of  R^.  From  our  point  of  view,  all  that  matters  is  the  finite-dimensional  space  of 
parameters  associated  with  the  positions  and  orientations  of  the  R-,  P-  and  S-joints.  Each  body 
serves  to  fix  the  relative  geometry  of  the  joints  attached  to  it.  Using  this  criterion  as  an 
equivalence  relation,  we  turn  attention  to  equivalence  classes  of  families.  We  observe  that  the 
equivalence  class  of  a  complete  family  is  a  finite-dimensional  space.  From  now  on  we  use  the 
term  "family"  to  denote  the  term  "equivalence  class  of  famihes,  the  latter  being  a  more  pre¬ 
cise  description  of  our  object  of  attention. 

There  is  an  important  special  class  of  families  that  merit  special  attention. 

Definition:  A  family  is  called  complete  if  all  the  parameters  of  the  above  four  types  are 
allowed  to  vary  independently. 

We  have  been  concentrating  on  the  theory  of  spatial  mechanisms,  but  all  of  our  statements  so 
far  apply  equally  well  to  planar  mechanisms.  The  adjustment  required  is  to  leave  out  spheri¬ 
cal  joints  and  consider  only  bodies  in  the  plane  z=0  connected  by  R-  and  P-  joints  that 
preserve  that  plane.  Now  that  we  are  considering  families  of  mechanisms,  however,  we  must 
recognize  the  distinction  between  the  planar  and  spatial  cases. 

Example  la:  Consider  a  four-bar  linkage  in  the  plane.  There  are  four  real  parameters:  the  dis¬ 
tances  between  the  revolute  joints  on  each  of  the  four  bodies.  In  this  case,  the  complete 
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family  is  the  fouT-dimensional  family  of  all  four-bar  linkages  in  the  plane. 


Example  lb:  Consider  a  three-bar  spatial  mechanism  with  spherical  joints  connecting  each  bar 
with  the  other  two.  There  are  three  real  parameters:  the  distances  between  the  S-joints  on  each 
of  the  three  bodies.  In  this  case,  the  complete  family  is  the  three-dimensional  family  of  all 
spatial  mechanisms  of  this  type. 

Note  that  the  four-dimensional  family  in  example  la  is  not  a  complete  family  when  the  4-bar 
is  viewed  as  a  spatial  mechanism.  In  fact,  unless  the  geometric  parameters  of  the  mechanism 
are  perturbed  in  special  ways,  there  are  no  solutions  to  the  system  of  algebraic  equations 
describing  the  perturbed  mechanism.  The  angles  and  positions  of  the  pairs  of  fixed  axes  asso¬ 
ciated  with  the  R-joints  on  the  four  bodies  must  satisfy  compatibility  conditions  if  the  bodies 
are  to  fit  together  in  any  way  at  all.  We  will  return  to  this  discussion  later. 

There  are  practical  reasons  for  considering  families  of  mechanisms  rather  than  individual 
mechanisms,  though  we  do  not  elaborate  this  point  here.  From  a  theoretical  point  of  view  we 
shall  see  that  the  "family  perspective"  helps  to  clarify  some  of  the  peculiar  mathematical 
features  of  special  mechanism  types. 

Before  proceeding,  we  require  some  more  definitions.  We  consider  the  planar  and  spatial 
cases  in  parallel. 


For  a  joint  j,  we  define  the  number  of  degrees  of  freedom  Dj  of  that  joint.  In  both  the  planar 
and  spatial  cases,  both  the  R-joint  and  the  P-joint  have  a  single  degree  of  freedom.  In  the 
spatial  case  only,  an  S-joint  has  three  degrees  of  freedom. 


In  the  following  we  let  Nj  denote  the  number  of  constrains  imposed  by  joint  j. 


10 


In  the  planar  case  Nj  =  3  -  Dj.  Each  body  in  the  plane  starts  with  three  degrees  of  freedom, 
but  this  number  is  reduced  by  Nj  independent  algebraic  equations  relating  the  positions  and 
orientations  of  the  two  bodies  connected  by  j. 

In  the  spatial  case  Nj  =  6  -  Dj.  Each  body  in  space  starts  with  six  degrees  of  freedom,  but 
this  number  is  reduced  by  Nj  independent  algebraic  equations  relating  the  positions  and  orien¬ 
tations  of  the  two  bodies  connected  by  j. 

A  little  thought  leads  to  the  following  formula  for  DOF,  where  DOF  is  defined  to  be  the 
dimension  of  the  configuration  space  of  a  mechanism: 

Planar  Mechanism:  DOF  =  3  (N-1)  -  S  Nj  (2.1a) 

j  e  J 


Spatial  Mechanism: 


DOF  =  6  (N-1)  -  S  Nj 
jeJ 


(2.1b) 


The  number  (N-1)  appears  because  one  of  the  bodies  is  fixed  to  ground. 

A  little  more  thought  reveals  that  these  formulas  do  not  always  work.  Considering  example  la 
as  a  spatial  mechanism  provides  one  counterexample  but  it  is  of  a  very  special  type.  We  pro¬ 
vide  another,  more  general  counterexample  to  clarify  the  problem. 

Example  2:  Consider  a  spatial  mechanism  consisting  of  two  bodies,  each  with  two  ball  joints. 
The  two  bodies  are  constructed  so  that  that  their  ball  joints  are  separated  by  the  same  dis¬ 
tance,  so  the  mechanism  can  be  assembled  in  R^.  Note  that  DOF  =  1,  because  the  ungrounded 
body  is  free  to  spin  about  the  axis  through  the  two  ball  joints.  The  right-hand  side  of 
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equation  (2.1b)  is  0. 


In  the  following  we  restrict  our  discussion  again  to  spatial  mechanisms.  The  reader  should  be 
able  to  make  the  proper  adjustments  for  the  restricted  planar  case. 

The  formula  of  equation  (2.1b)  does  not  work  for  example  2  because  the  mechanism  type  is 
overconstrained.  By  overconstrained,  we  mean  that  a  special  relation  must  exist  among  the 
parameters  of  the  complete  family  in  order  for  the  mechanism  to  be  realizable.  If  the  distances 
between  ball  joints  on  the  two  bodies  are  different,  no  solutions  (real  or  complex)  exist  to  the 
system  of  algebraic  equations  describing  the  mechanism. 

To  address  this  situation,  we  consider  the  configuration  space  C  of  the  mechanism  as  a  subset 
of  the  product  space  E(3)^"^^  cut  out  by  the  intersection  of  the  constraint  equations  [recall 
that  one  of  the  bodies  is  grounded,  hence  the  exponent  (N-1)]. 

Definition:  A  mechamsm  is  called  simply-constrained  if  its  configuration  space  C  is  the 

transverse  intersection  of  the  L  N;  algebraic  constraint  equations  imposed  by  its  joints. 

jeJ 

For  a  spatial  mechanism,  the  transversality  condition  means  that  the  rank  of  the  subspace  of 

the  cotangent  space  of  E(3)^"^^  spanned  by  the  differentials  of  the  constraint  equations  at 

each  point  of  C  is  equal  to  E  Nj.  From  this  definition  the  following  consequences  are 

j  e  J 

immediate  for  a  simply-constrained  mechanism: 

1)  The  configuration  space  C  is  a  smooth  maitifold 

2)  The  dimension  of  C  is  DOF,  which  is  6  (N-1)  -  E  Nj 

3)  A  small  neighborhood  of  the  complete  family  is  realizable  by  simply  constrained  mechanisms 
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Now  consider  the  parameter  space  X  of  the  complete  family  of  a  simply-constrained  mechan¬ 
ism  M.  The  mechanism  M  represents  a  point  in  X  at  which  the  constraint  equations  intersect 
transversely  in  This  transversality  condition  is  an  open  condition,  hence  is  true  for  a 

Zariski-open  set  U  in  X.  The  mechanisms  corresponding  to  points  in  U  are  simply- 
constrained.  Therefore,  it  makes  sense  to  speak  of  a  family  of  simply-constrained  mechan¬ 
isms;  we  call  such  a  family  a  simple  family.  Note  that  not  every  mechanism  in  a  simple  fam¬ 
ily  need  be  simply-constrained  -  in  fact,  there  will  be  closed  sets  in  X  for  which  the 
configuration  spaces  are  singular  algebraic  sets.  On  the  other  hand,  every  Zariski-open  set  of 
mechanisms  containing  one  of  these  singular  algebraic  sets  will  contain  simply-constrained 
mechanisms,  so  the  singular  behavior  inside  a  simple  family  is,  in  some  sense,  not  so  bad. 

The  difference  between  example  lb  and  2  is  now  clear  ~  example  lb  lies  in  a  simple  family 
while  example  2  does  not.  A  mechanism  that  does  not  lie  in  a  simply  family  is  called  over¬ 
constrained. 

The  general  concept  of  mechanism  families  seems  natural  from  a  mathematical  viewpoint, 
though  our  current  understanding  of  them  is  quite  limited.  In  his  undergraduate  thesis 
[Walker],  Walker  analyzed  the  class  of  simple  families  of  single-loop  planar  linkages.  For 
this  class  the  complete  families,  when  normalized  by  a  geometric  scale  factor,  are  geometric 
simplexes  (higher-dimensional  tetrahedra)  of  dimension  one  less  than  the  number  of  bars  in 
the  linkage.  Walker  determined  the  structure  of  the  algebraic  subsets  corresponding  to  singu¬ 
lar  mechanisms  within  the  families,  and  constructed  Morse  functions  on  the  corresponding 
real  configuration  spaces  to  determine  the  changes  in  topology  for  fanulies  of  mechanisms 
passing  through  those  singular  sets.  He  also  completed  the  classification  of  two-DOF 
configuration  spaces  for  mechanisms  of  this  type,  and  obtained  partial  results  for  the  higher¬ 
dimensional  cases.  The  local  parametrization  problem  appears  tractable  for  planar  linkages  of 
arbitrary  loop-structure,  though  Walker  did  not  address  it.  We  are  not  aware  of  any  other 
references  along  these  lines. 


13 


2.4  Statement  and  Discussion  of  the  Main  Problem 

We  are  interested  in  the  following  problem: 

MAIN  PROBLEM:  Determine,  by  algebraic  methods,  explicit  parametrizations  of  the 
configuration  space  for  a  family  of  specified  mechanisms.  This  parametrization  should  be  of 
minimal  degree. 

The  rest  of  this  section  is  devoted  to  discussion  of  the  main  problem. 

As  we  have  already  seen,  for  most  mechanisms  in  a  simple  family  (i.e.  a  Zariski-open  set), 
the  configuration  spaces  are  smooth  manifolds  of  dimension  DOF  as  defined  by  equation  (lb). 
When  C  is  a  manifold,  the  explicit  parametrizations  we  have  in  mind  form  a  special,  algebraic 
coordinate  atlas  for  C. 

In  general  (simply-constrained  or  not),  the  construction  is  as  follows. 

From  the  set  J  of  all  joints,  select  (if  possible)  a  subset  Q  such  that: 

1)  Z  D:  =  DOF 
jeQ 

2)  the  degrees  of  freedom  of  joints  in  Q  are  independent 

For  many  practical  examples  we  have  analyzed,  subsets  Q  satisfying  the  first  condition  do 
exist  (see  example  2  above  for  a  counterexample).  The  meaning  of  the  second  condition  will 
become  apparent  immediately. 

Consider  the  situation  where  all  the  free  parameters  associated  with  the  joints  in  Q  are  frozen 
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at  some  fixed  values.  When  C  is  a  smooth  manifold,  the  two  conditions  above  are  those 
required  to  make  the  corresponding  set  of  configurations  a  discrete,  zero-dimensional  set  (i.e. 
a  set  of  points).  This  discrete  set  is  the  zero  set  of  the  ideal  generated  by  finitely  many  poly¬ 
nomial  equations  and  so  is  a  finite  set.  The  number  of  points  in  that  set  (counting  multiplici¬ 
ties)  is  denoted  Dcvq,  the  degree  of  C  over  Q.  It  is  not  always  the  case  that  sets  Q  can  be 
found  satisfying  the  above  two  conditions  even  for  simple  families  (e.g.  example  lb),  but 
there  are  enough  examples  to  make  this  special  condition  worth  consideration. 

For  the  rest  of  this  subsection  we  assume  a  set  of  joints  Q  satsifying  the  above  two  conditions 
can  be  found. 

The  geometric  situation  can  viewed  in  another  way  as  follows.  Suppose  the  configuration 
space  C  of  the  mechanism  M  is  a  smooth  manifold  and  let  Q  be  a  set  of  joints  as  above.  Let 
Hq  be  the  product  manifold  of  all  the  geometric  subgroups  of  E(3)  corresponding  to  the  rela¬ 
tive  motions  of  each  joint  in  Q.  Then  there  is  a  finite  branched  covering  map  tc  :  C  — »  Hq 
[Gunning].  The  fiber  over  each  point  h  of  Hq  is  the  finite  set  of  points  in  C  for  which  the 
joints  in  Q  have  the  values  specified  by  h.  The  number  Dq/q  is  the  branching  order  of  this 
covering.  The  set  defined  to  be  the  union  of  multiple  points  in  C  with  respect  to  it  is 
called  the  branch  locus  of  tc. 

Now  we  observe  that  the  construction  above  for  the  configuration  space  C  of  a  single  mechan¬ 
ism  extends  naturally  to  the  configuration  spaces  of  all  mechanism  in  a  family  (assuming  the 
family  is  realizable,  as  in  the  case  of  simple  families).  For  any  configuration-space  manifold 
C  in  that  family  the  number  Dcvq  is  the  same.  Because  the  number  Dc'/q  is  independent  of 
the  choice  of  O',  we  change  notation  and  denote  this  number  by  Dq. 

We  are  finally  able  to  provide  a  precise  statement  of  the  main  problem.  Consider  a  family  of 

mechanisms  and  a  subset  Q  satisfying  the  above  conditions  1  and  2.  Determme  a  mapping  (in 
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the  form  of  a  computational  algorithm)  from  the  product  space  Hq  to  a  set  of  polynomials  in 
flag  form  for  the  remaining  joint  positions  [Cox],  That  flag  form  should  have  the  property 
that,  for  a  Zariski-open  set  of  mechanisms  in  the  family,  and  for  every  fixed  point  in  Hq,  the 
number  of  discrete  solutions  (counting  multiplicity)  in  that  flag  is  Dq, 

The  motivation  for  the  main  problem  arises  in  the  fi-amework  of  computer-aided  mechanism 
design,  analysis  and  djmamic  simulation.  From  the  flag  form  it  is  possible  to  compute 
efficiently  all  allowed  positions  of  a  mechanism.  By  constructing  the  flag  for  a  family  of 
mechanisms,  we  can  use  the  same  basic  algorithm  for  mechanisms  in  a  general  class,  thereby 
allowing  a  mechanism  designer  to  change  the  dimensions  (e.g.  lengths)  of  the  bodies,  the 
orientations  and  positions  of  the  joints,  in  a  interactive  manner.  For  simple  families,  small 
changes  can  be  made  to  every  one  of  the  parameters  while  staying  within  the  class  of  simply- 
constrained  mechanisms.  Our  computational  procedure  provides  a  good  coordinate  atlas  for 
all  mechanisms  in  that  local  family. 

Once  the  algebraic  equations  are  formulated,  a  polynomial  flag  of  the  type  just  described  can 
be  computed  by  the  Groebner-Basis  algorithm  [Buchbergerl,2].  In  practice,  however,  the 
Groebner-Basis  algorithm  does  not  provide  answers  quickly  (see  [EM]  for  a  discussion  of 
computational  issues).  The  subject  of  this  paper  is  an  approach  that  produces  a  flag  of  the 
same  type  by  more  practical  computational  methods. 

In  general,  a  single  subset  Q  will  not  provide  a  good  atlas  for  an  entire  configuration  manifold 
because  the  projection  map  jc  generally  has  a  non-empty  branch  locus.  To  form  a  good  atlas, 
one  may  consider  a  finite  collection  of  subsets  Qi,  •  •  •  ,Qk  for  which  the  intersection  of  the 
corresponding  branch  loci  is  empty.  A  good  atlas  is  needed  for  working  on  dynamical  prob¬ 
lems  such  as  mechanism  control. 
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3.  DEVELOPMENT  OF  THE  MULTI-AFFINE  SYSTEM 


For  a  given  mechanism,  each  of  the  rigid  bodies  in  W  has  its  own  reference  coordinate  sys¬ 
tem.  From  a  kinematic  viewpoint,  the  role  of  the  body  coordinate  system  is  to  allow 
specification  of  the  geometric  relation  (separation,  relative  orientation  in  space)  among  the 
joints  on  that  body.  The  mathematical  data  describing  the  body  in  its  coordinate  chart  is  a 
subset  of  the  data  required  by  a  machinist  to  build  it.  Within  each  coordinate  chart  a  three- 
dimensional  coordinate  system  is  specified,  and  the  coordinates  of  salient  parts  of  the  body  are 
identified.  In  the  following  we  will  use  the  notation  Uj  to  denote  the  coordinate  chart  of  body 
j,  and  (Xj,  yj,  Zj)  to  denote  the  coordinate  functions  in  Uj.  When  discussing  properties  that  hold 
for  a  general  chart  we  often  omit  the  subscript  j. 

3.1  Loop  Closure  Equations 


We  first  discuss  the  notion  of  transition  functions.  Roughly  speaking,  a  transition  function  is  a 
mapping  from  one  coordinate  chart  to  another  that  identifies  those  points  that  are  joined 
together  in  the  assembled  mechanism.  In  the  next  paragraph  we  clarify  this  notion. 

Let  SO(3)  denote  the  special  orthogonal  group  of  real  orthogonal  matrices  of  determinant  1. 
The  Euclidean  group  E(3)  is  the  group  of  all  mappings  ^  from  to  R^  of  the  form: 

<l)(x)  =  A  X  +  b  (3.1) 

where  A  is  in  SO(3)  and  b  is  in  R^.  There  is  a  standard  representation  of  E(3)  as  a  subgroup 
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of  4  X  4  matrices  by  the  mapping  p  defined  as  follows: 


b 

1 


(3.2) 


where  <(),  A  and  b  are  as  in  equation  1.  The  value  of  the  transformation  <1>  applied  to  the 
three- vector  x  is  realized  by  identifying  x  with  the  four- vector  x',  where 


X  = 


X 

1 


(3.3) 


and  computing  the  usual  matrix  product 

[<t)(x)]'  =  p(<|))  x' 


(3.4) 


Consider  two  coordinate  charts,  Uj  and  U2,  for  a  pair  of  bodies  connected  by  a  joint.  Using 
the  4  X  4-matrix  notation,  the  set  of  Euclidean  motions  <>12  satisfying  the  constraint  equations 
have  the  following  forms,  depending  on  joint  type. 

First,  for  an  R-type  joint,  it  is  easy  to  check  that  the  1-parameter  transition  functions  are  of 
the  form 


P(<|)2i(e))  = 


R(e)  b-R(0)a 
0  1 


(3.5) 


where  a,  b  are  constant  vectors  that  represent  fixed  points  of  the  relative  motion  between  the 
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two  bodies  in  the  two  coordinate  systems.  The  matrix  R(0)  is  a  rotation  matrix  of  the  form: 


(3.6) 


for  0  an  arbitrary  angle.  V  the  axis  of  rotation  in  U2,  W  the  axis  of  rotation  in  Ui,  Vp^rp  an 
orthonormal  2x3  complement  to  V,  and  an  orthonormal  2x3  complement  to  W. 


For  a  P-type  joint,  the  one-parameter  transition  functions  are  of  the  form: 


P(<|)2i(t))  = 


I3 

0 


t  a 
1 


(3.7) 


where  t  is  a  scalar  variable  and  a  is  a  constant  vector  in  the  direction  of  the  motion.  The 
matrix  I3  is  the  3x3  identity  matrix. 

For  an  S-type  joint,  the  most  general  transition  function  is  of  the  form: 


P«1)2i(R))  = 


R  b-Ra 
0  1 


(3.8) 


where  a,  b  are  constant  vectors  that  represent  fixed  points  of  the  relative  motion  between  the 
two  bodies  in  the  two  coordinate  systems.  R  is  a  general  rotation  matrix. 

Now  that  the  transition  functions  have  been  defined  it  is  an  easy  matter  to  write  down  the 
loop-closure  equations  (LCEs).  The  loop-closure  equations  are  obtained  by  taking  the 
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composition  of  transition  functions  around  a  closed  loop.  Starting,  say,  in  coordinate  chart  Ui, 
pick  a  general  point  Xj.  Any  transition  function  <t)2i  maps  to  its  image  X2  in  the  coordinate 
chart  U2  (change  of  coordinates).  The  point  X2  is  mapped  in  turn  to  X3  in  U3  by  a  transition 
function  <1)32,  and  so  on,  until  finally  the  point  Xl  in  Ul  is  mapped  by  the  function  <t)iL  to  y^  in 
the  original  chart  Uj.  When  the  parameters  of  the  transition  functions  are  those  one  would 
measure  on  the  mechanism  in  a  realizable  configuration,  the  compatibility  relation  yj  =  xj  is 
satisfied.  This  compatibility  relation  holds  for  all  xj  in  Ui  for  which  the  composition  is 
defined.  This  relation  is  represented  mathematically  by  a  loop  closure  equation. 

For  a  realizable  configuration,  the  con^atibility  relation  must  hold  for  the  following  geometric 
reason.  Viewing  the  transition  functions  as  changes  of  coordinates  associated  with  the 
different  Ixxlies,  we  see  that  yx  and  xx  are  the  coordinates  of  the  same  geometric  point  in  the 
same  coordinate  chart  Ux,  so  they  are  equal.  The  associated  loop  closure  equation  states  that 
the  composition  is  the  identity  map  in  E(3). 

Consider  the  example  of  the  7R  mechanism.  Because  this  mechanism  has  only  a  single  loop, 
all  l<x)p  equations  have  the  form: 

<l)j+7  j+6(ej^.6)  <t)j+6  j+sCOj+s)  •  •  •  <t>j+i  j(ej)  =  identity  (3.9) 

where  the  subscripts  are  evaluated  mod  7.  As  an  example,  one  of  these  equations  in  terms  of 
the  4  X  4  matrices  is  (pick  coordinates  so  that  bj  =  0): 

^17(67)  "^17(07)  ^7 
0  1 

where  I  is  the  4  x  4  identity  matrix,  Rj+x  j(0j)  is  the  subscripted  version  of  equation  6,  and  aj 


R32(®2)  "^32(^2)  ^2 
0  1 


R2i(0i) 

0 


”R2i(®i) 

1 


=  I  (3.10) 
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is  the  subscripted  version  of  equation  (3.5). 

The  readCT  familiar  with  differential  geometry  will  recognize  the  significance  of  the  above 
construction  in  more  concise  but  abstract  terms.  What  we  have  described  is  a  special  type  of 
flat-Euclidean  manifold  structure  [Charlap]  on  a  geometric  space  associated  with  the  mechan¬ 
ism.  This  abstract  viewpoint  was  the  original  motivation  behind  the  transition  function 
approach.  This  abstract  point  of  view  does  not  tell  us  how  to  proceed  with  the  computations, 
but  it  did  get  us  started.  The  quaternion  pair  notation  described  next  provides  the  computa¬ 
tional  key. 

3.2  The  Quaternion  Pair  Notation 


Quaternions  are  often  used  in  kinematic  analysis,  so  the  mechanism  community  is  familiar 
with  their  properties.  Even  so,  it  appears  that  the  quaternion-based  approach  we  present  here 
has  not  been  exploited  in  mechanism  analysis. 

The  point  of  introducing  quaternions  is  to  transform  the  loop  closure  equations  (LCEs) 
derived  in  Section  3.1  into  a  simpler  form.  The  simplification  is  not  merely  cosmetic,  it  is  a 
basic  reduction  of  the  system  of  equations  into  an  equivalent  but  more  tractable  form. 

We  assume  the  reader  is  familiar  with  the  basics  of  quaternions.  A  general  reference  is  [Por- 
teus].  The  section  below  is  a  slight  generalization,  of  the  more  detailed  presentation  in  [ME], 
so  some  details  are  omitted. 

3.2.1  Unit  Quaternions  and  SO(3) 


Tha-e  is  a  well-known  covering  map  from  the  unit  quaternions  Q(l)  to  the  group  of  rotation 
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matrices  S0(3).  This  brief  section  reviews  that  construction  and  some  of  its  immediate  corol¬ 
laries. 

Let  p  be  a  nonzero  quaternion.  Then  pwp*  is  a  pure  quaternion  if  and  only  if  w  is  a  pure 
quaternion.  By  the  mapping  w  pwp*  each  nonzero  quaternion  p  defines  an  invertible  linear 
transformation  S(p)  of  T.  We  will  primarily  consider  this  transformation  in  the  special  case 
where  p  is  a  unit  quaternion,  in  which  case  S(p)  is  in  SO(3).  Relative  to  the  basis  {i,j3c}  of  T 
(which  we  tacitly  assume  in  the  sequel),  the  matrix  S(p)  associated  with  p  has  the  form. 

Po+Pf  P1P2-P0P3  P1P3  +  P0P2 
S(Po.  Pi.  P2.  P3)  =  2  P1P2  +  P0P3  Po  +  pI  P2P3  ~  PoPi 
P1P3-P0P2  P2P3  +  P0P1  P0+P3  . 

The  mapping  p  S(p)  projects  the  unit  quaternions  Q(l)  onto  the  rotation  matrices  SO(3). 
For  each  matrix  K  in  SO(3)  there  are  two  quaternions,  p  and  -p,  such  that  S(p)  =  S(-p)  =  K. 
Furthermore,  S(p2Pi)  =  S(p2)S(pi). 

The  geometric  link  between  rotation  matrices  and  their  associated  unit  quaternions  is  easily 
derived.  Specifically,  any  unit  quaternion  q  has  the  form: 

q  =  cos(y)  +  sin(Y)  w  (3.12) 

where  w  is  a  unit-length,  pure  quaternion.  The  matrix  S(q)  is  a  rotation  by  the  angle  2y  about 
the  axis  in  the  w-direction. 

3.2.2  The  Quatenion-Pair  Cover  of  the  Euclidean  Group 


1  0  0 

-010.  (3.11) 

P  0  1. 
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Now  we  turn  attention  to  a  special  group  defined  by  quaternion  pairs.  Specifically,  we  will 
look  at  the  set  Q(l)  x  T  of  ordered  pairs  of  quaternions  (q,w)  and  define  a  group  structure  on 
this  set.  Let  qi,  q2  be  unit  quaternions  and  Wj,  W2  be  pure  quaternions.  Define  the  product: 

(q2.  W2)(qi»  wi)  =  (q2qi.  W2  +  q2Wiq2*)  (3.13) 

Observe  that  the  product  is  well-defined  in  Q(l)  x  T.  Note  that  (1,0)  is  the  identity  element, 
and 


(q,  w)  ^  =  (q*,  -q*wq) 


(3.14) 


The  set  Q(l)  x  T  with  this  product  forms  a  group,  which  we  call  the  quaternion-pair  group. 
This  group  is  the  universal  cover  of  the  Euclidean  group  E(3).  The  full  power  of  the  covering 
space  stracture  is  not  needed  for  our  application,  all  we  need  is  the  definition  of  the  mapping. 
Therefore,  we  construct  the  covering  map  S*  below,  but  we  do  not  address  its  topological  pro¬ 
perties. 

The  covering  map  from  the  quaternion-pair  group  to  E(3)  is  closely  related  to  the  map  S 
defined  in  Section  3.2.1  that  maps  Q(l)  onto  SO(3).  Let  <|)i,  02  be  two  elements  of  E(3)  and 
suppose  for  j  =  1,2: 


P(0j)  = 


0  1 


(3.15) 


Pick  qj  such  that  S(qj)  =  Rj,  and  let  wj  be  the  pure  quaternion 
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(3.16) 


Wj  =  aj.ii  +  ^21  +  aj,3k 


where  aj  =  (aj,!.  a^,  aj.3).  Defining  the  map  S*[(qj,  Wj)]  =  <|)j,  it  is  an  easy  computation  to  ver¬ 
ify  that 


S*[(q2.  W2)(qi,  Wi)]  =  ^2  <t>i 

The  mapping  S*  is  clearly  onto,  it  is  a  two-to-one  mapping  like  k. 


3.3  Transition  Functions  in  Quaternion-Pair  Form 


The  primary  goal  of  this  subsection  is  to  determine  a  quaternion-pair  form  for  transition  func¬ 
tions  equivalent  to  the  expressions  in  equations  (3.5)  through  (3.8)  in  Section  3.1.  The  idea  is 
to  represent  the  transition  function  (we  ignore  the  dependence  of  (()  on  its  parameters) 
by  a  quaternion  pair: 


j)  =  (qj»  Wj)  (3-18) 

where  qj  is  a  unit  quaternion  and  Wj  is  a  pure  quaternion.  We  begin  with  the  special  case  of 
R-joints. 

3.3.1  Representing  R-joints 


Consider  the  expression  in  equation  (3.5),  and  let  <|)  denote  <(>21.  Fix  0  and  select  q,  one  of  the 
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two  quaternions  such  that  S(q)  =  R(0).  Let  a  ,  P  be  the  pure  quaternions  associated  with  the 
vectors  a  and  b.  The  results  of  the  previous  subsection  suggest  that  we  use 

=  (q.  P  -  qaq*)  (3*19) 


for  the  quaternion  pair  representing  (|).  Furthermore,  if  <t)i  and  ^2  ^  two  transition  functions 
corresponding  to  (qi,  pj  -  qittiqj*)  and  (qj,  P2  -  q2a2q2*)*  then  the  composition  map  is  the 
product  (see  equation  (3.19)): 

o(<l)2<l>i)  =  (qiqi.  p2  +  q2(Pi  -  -  q2qi«iqi*q2*)  (3-20) 

There  is  a  pattern  in  expression  on  the  right  hand  side  of  this  last  equation  that  will  be 
clarified  in  Section  3.4. 

A  vital  detail  in  the  quatemion-pair  representation  is  to  include  the  0-dependence,  which  is 
described  in  the  following  circle  lemma. 

Circle  Lemma:  Let  R(0)  be  the  one-parameter  family  of  rotation  matrices  defined  in  equation 
(3.6).  There  is  a  pair  of  unit  quaternions  (p,q),  with  p  orthogonal  to  q,  such  that  the  great  cir¬ 
cle  Cj  in  Q(l)  defined  by 


Cl  =  {cos(y)  P  +  sin(Y)  q  I  0^  Y  -  27t} 


(3.21) 


is  mapped  by  S  onto  the  circle  C2  in  SO(3)  defined  by 


C2  =  {R(0)  I  0^  0  ^  271} 


(3.22) 
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The  proof  of  this  lemma  is  found  in  [ME]. 


A  more  geometric  statement  of  the  circle  lemma  is:  circles  in  the  rotation  group  are  covered 
by  great  circles  in  the  unit  quaternion  group.  The  significance  of  this  result  will  be  realized  in 
Section  3.4. 

3.3.2  Representing  S-joints 

Consider  the  expression  in  equation  (3.8),  and  let  (j)  denote  Select  q,  one  of  the  two 
quaternions  such  that  S(q)  =  R  and  let  a  ,  p  be  the  pure  quaternions  associated  with  the  vec¬ 
tors  a  and  b.  Then  as  above  we  choose: 

Cf(<l))  =  (q,  P  -  qaq*)  (3.23) 

for  the  quaternion  pair  representing  <t).  The  formulas  for  products  of  transition  functions  asso¬ 
ciated  with  both  S-joints  and  R-joints  are  the  same  as  equation  (3.20)  above.  The  only 
difference  is  that  the  variable  q  is  now  an  arbitrary  unit  quaternion,  and  not  of  the  special 
form  stated  in  the  circle  lemma. 

3.3.3  Representing  P-joints 

This  is  the  easiest  case.  For  <|)(t)  as  in  Equation  (3.7)  the  answer  is: 

a((l)(t))  =  (qo,t  a)  (3.24) 

where  qo  is  a  constant  unit  quaternion  and  a  is  the  pure  quaternion  associated  with  the  vector 
a. 
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3.4  The  Multi-Affine  Equations 


Examining  the  quatemion-pair  representations  described  in  the  last  section  we  find  that  the 
LCEs  have  the  following  form: 


(qi,  wi)  •  •  •  (qL,  Wl)  =  (1,0)  (3.25) 

We  now  observe  a  particular  feature  of  the  expression  on  the  LHS  of  this  last  equation: 

1)  the  unit-quaternion  part  of  a  product  is  the  product  of  the  unit  quaternions 

2)  each  summand  of  the  pure  quaternion  part  is  a  sequential  conjugate  expression 

By  sequential  conjugate  expression,  we  mean  an  expression  of  the  form: 

F(w)  =  Qi  •  "  Qm  w  qn,*  •  •  •  Qi*  (3.25) 

where  w  is  a  pure  quaternion  and  m  varies  from  0  to  L.  The  point  is  that,  for  all  m,  when  a 
sequential  conjugate  expression  is  right-multiplied  by  the  unit-quatemion  part  of  the  product, 
the  result  is  multilinear  in  all  the  unit  quaternions  qj. 

In  addition,  the  equation: 


Qu(qi  •  •  •  qj  =  0 


(3.27) 


where  the  function  Qu  means  "pure-quatemion  part,"  provides  another  set  of  multi-linear 
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equations  in  the  variables  qj. 


These  multi-linear  equations  in  qj  can  be  reduced  to  multi-affine  equations  in  fewer  variables 
by  the  following  devices: 

1)  For  each  R-joint  in  the  product,  use  the  Circle  Lemma  and  divide  by  cos  (y) 

2)  For  each  S-joint  in  the  product,  divide  by  qg  =  Re(q) 

0i 

Let  Zj  denote  the  variable  tan(-— ),  and  define  the  variables  ej,  fj,  gj  by: 


ej 

_  j[_ 

qj.i 

h 

Sj_ 

qj.o 

%2 

(3.28) 


Then  the  LCEs  become  multi-affine  in  the  variables  zj,  ej,  fj,  gj. 

Observe  that  the  LCEs  are  already  multi-affine  in  the  variables  ^  associated  with  P-joints.  The 
discussion  just  concluded  brings  us  to  the 

Main  Result:  The  loop  closure  equations  (LCEs)  generate  a  host  of  multi-affine  equations, 
forming  a  modular  system  of  a  very  special  type.  The  number  of  independent  variables 
Zi,  Cj,  fj,  gj  tk  appearing  in  these  equations  is  equal  to  the  sum  over  all  joints  of  the  number 
of  joint  degrees  of  freedom.  In  general,  there  are  far  more  equations  than  unknowns. 

The  primary  benefit  of  this  result  lies  in  providing  a  special  structure  to  the  set  of  equations 
whose  solutions  describe  the  configuration  spaces  of  mechaiusms  of  this  type.  An  additional 
feature  is  the  systematic  approach  it  provides  for  generating  these  equations. 
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An  approach  for  numerical  solution  of  multi-affine  equations  of  this  type  is  illustrated  in  the 
solution  to  the  7R  problem  presented  in  [ME]. 
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4.  SUMMARY 


The  discussion  below  is  a  summary  of  our  six-step,  general  approach  for  mechanism 
kinematic  analysis.  The  details  for  a  particular  application  are  presented  in  [ME]. 

We  have  seen  that  the  quaternion-pair  approach  described  below  can  be  applied  to  spherical, 
prismatic  and  cylindrical  joints  as  well  as  the  revolute  joints  illustrated  in  [ME].  Systems  of 
multi-affine  equations  are  obtained  in  the  same  way.  Multi-loop  mechanisms  can  also  be  han¬ 
dled.  One  of  the  primary  goals  of  our  research  effort  to  use  the  theory  on  problems  of 
moderate  size  to  determine  the  practical  (computational)  limitations  of  this  genCTal  approach 
on  more  complex  mechanisms. 

4.1  Reference  Coordinate  Systems 

Each  rigid  piece  of  the  mechanism  has  its  own  reference  coordinate  system.  The  concept  is 
simple  -  each  piece  can  be  thought  of  as  a  separate  entity  standing  by  itself  in  Euclidean 
three-space. 

4.2  Transition  Functions 


The  transition  functions  represent  the  translation  and  rotation  (Euclidean  motion)  associated 
with  the  change  of  reference  coordinate  systems  from  one  body  to  the  next.  The  transition 
function  contains  the  fixed  parameters  of  the  body  as  weU  as  the  free  variables  associated  with 
the  joint. 
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4.3  Quatemion-Pair  Notation 


The  quatemion-pair  notation  consists  of  a  pair  of  quaternions  that  together  represent  a  single 
Euclidean  motion;  the  first  is  a  unit  quaternion  that  represents  the  rotation  and  the  second  is  a 
pure  quaternion  that  represents  translation.  The  advantage  of  this  notation  is  that  it  leads  to 
loop  equations  that  are  multi-affine  in  the  joint  variables.  These  multi-affine  equations  are  of 
lower  total  degree  and  therefore  are  much  easi^  to  solve  than  the  multi  quadratic  equations 
that  arise  from  using  the  standard  4x4  matrices  containing  rotation  matrices  and  a  translation 
vector. 

The  quaternion-pairs  form  a  mathematical  group  under  a  product  defined  in  Section  3.2.2.  The 
composition  of  transition  functions  is  obtained  by  multiplying  the  corresponding  quaternion 
pairs. 

4.4  Loop  Equations 

The  loop  closure  equations  are  obtained  by  setting  the  product  of  all  the  quatemion-pair  tran¬ 
sition  functions  around  a  loop  equal  to  the  quatemion-pair  identity. 

4.5  Multi-Affine  Equations 

After  multiplication  by  the  appropriate  factors,  the  resulting  loop  equations  are  multi-affine  in 
the  joint  variables.  Because  the  quatemion-pair  multiplication  is  not  commutative,  additional 
independent  equations  can  be  formed  by  cyclic  permutations  of  the  transition  functions  (start¬ 
ing  at  a  different  body  when  foiming  the  loop).  Some  of  the  equations  obtained  in  this  way 
involve  fewer  variables,  so  more  multi-affine  equations  of  the  same  degree  can  be  formed  by 
multiplying  those  equations  by  each  of  the  missing  variables.  In  this  way,  a  large  number  of 
affine  equations  can  be  formed. 
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4.6  Elimination  Techniques 

Standard  polynomial  elimination  techniques  (such  as  Groebner  basis)  work  on  general  systems 
of  polynomials,  but  their  run  times  are  doubly  exponential  in  the  number  of  noninteger  param¬ 
eters.  This  makes  them  unusable  for  mechanisms  of  even  moderate  complexity.  Conse¬ 
quently  we  developed  specific  elimination  techniques  based  on  matrix  polynomial  methods 
[Wedderbum]  that  are  much  faster  for  the  specific  (multi-affine)  type  of  systems  of  polynomi¬ 
als  arising  from  mechanism  kinematics.  These  specific  techniques  consist  of  generalized 
eigenvalue  problems,  explicit  polynomial  representation  of  the  null  space  of  affine  matrices, 
and  numerical  methods  fi'om  linear  algebra. 
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Abstract 


This  paper  deals  with  n  =  wi +712 H" •••■!* equations,  multi-affine  in  k  sets 
of  variables  with  n»  variables  in  the  i*'*  set.  In  [St]  a  formula  is  given  for,  d, 
the  number  of  solutions,  and  a  degree  d  polynomial  in  one  variable  is  formed 
from  a  hyper-determinant  function  of  the  coefficients  and  one  variable.  The 
d  roots  of  this  polynomial  are  the  value  of  that  variable  at  the  d  solution 
points. 

In  the  current  paper,  we  give  an  alternative  proof  for  the  number  of  solu¬ 
tions,  and  we  form  a  d  x  d  eigenproblem  whose  eigenvalues  are  the  d  values  of 
one  variable,  and  the  other  variables  are  ratios  of  entries  of  the  eigenvectors. 
This  technique  of  forming  the  eigenproblem  directly,  then  solving  it  numer¬ 
ically  gives  an  0(cP)  algorithm  (where  d  »  n).  The  previous  technique 
[St], [Can],  [GKZ]  of  forming  the  degree  d  polynomial  is  0(d")  or  0((nA:)")  ? 

The  three  steps  in  our  algorithm  are: 

1)  Use  dilation  to  form  a  [mm(ni)d]  x  [min{ni)d],  rank-d,  eigenproblem. 

2)  If  min{ni)  >  1,  use  deflation  to  reduce  the  size  of  the  eigenproblem 

a)  Compute  svd  of  a  matrix  of  size  [min(ni)d]  x  d. 

b)  Compute  svd  of  a  matrix  of  size  [{min{ni)  —  l)d]  x  [min{ni)d\. 

3)  Solve  a  d  X  d  generalized  eigenproblem. 

The  matrices  involved  in  each  of  the  above  three  steps  are  either  linear 
in  the  coefficients  of  the  original  equations,  or  are  the  product  of  orthogonal 
matrices  times  matrices  that  are  linear  in  the  coefficients  of  the  original 
equations.  The  d  eigenvalues  give  the  d  values  of  one  of  the  variables,  while 
the  corresponding  values  of  each  of  the  other  variables  are  given  by  a  ratio  of 
two  linear  combinations  of  the  eigenvector.  The  runtime  of  the  algorithm  is 
dominated  by  the  runtime  of  the  svd’s  and  the  eigensolver.  The  total  number 
of  floating  point  operations  is  10[min{ni)d\^ . 

This  paper  also  deals  with  the  case  where  there  are  more  equations  than 
unknowns.  The  extra  equations  may  eliminate  some  or  all  of  the  solutions. 


Chapter  1 

Same  Number  of  Equations  as 
Unknowns 


1 . 1  Introduction 

Systems  of  multi-aflBne  equations  arise  from  applications  where  some  specified 
input  vector  Ui,  is  multiplied  by  a  sequence  of  k  subsystems  to  give  a  specified 
output  vector  For  i  =  1  to  A:,  subsystem  i  is  linear  in  its 

own  subset  of  rii  variables,  and  is  of  the  form: 

y,  =  Hi{^)ui  (1) 


where 


d"  •••  “I" 

(2) 

Let 

TUi+i  xmi  =  size(Hi{s^)) 

(3) 

To  ensure 

as  many  equations  as  unknowns  at  each  step: 

i 

Uj  for  i  <k 

(4) 

i=i 

and 

k 

“^k+l  =  X) 

(5) 

j=i 


1 


The  product  of  these  subsystems  then  gives: 


Ik 


i+i=Ui 


he = +  +  •••  +  2kl 


(6) 


(7) 


Given  the  value  of  the  output  vector  ^  G  i?("i+"2+  -+"»=)  and  the  value 
of  the  input  vector,  tti,  we  have  ni  +  fi2  +  ...  +  7ik  multi-affine  equations  in 
ni  +  712  +  •••  +  wjk  unknown  Xij. 


Applications  that  involve  systems  of  multi-affine  equations  include: 


Generalized  nxn  eigenvalue  problem,  where  A:  =  2,  ni  =  n  —  1,  n2  —  1) 
=  Q,  Hijui  =  column  j+l  of  the  n-dimensional  identity  matrix.  This 

gives: 


0  =  [^f2,0  +  2:2,1  i?2,l] 


a?i,i 

a;i,2 


(8) 


^l,n— 1 


Robotics  [ME96],  where  each  subsystem  is  a  link  in  a  robot  arm,  the 
variables  are  associated  with  the  joint  degrees  of  freedom,  the  inputs  are  the 
base  position  and  attitude,  and  the  outputs  are  the  end-effector  positions 
and  attitude. 

Electronic  circuits,  where  multi-input  multi-output  subsystem  i  contains 
Ui  op-amps,  the  variables  are  the  impedances  ratios  for  each  op-amp,  the 
input  is  a  vector  of  voltages  input  to  the  first  stage,  and  the  output  is  a 
vector  of  voltages  output  from  the  last  stage. 


In  ^sterns  like  these,  we  can  solve  for  all  values  of  each  parameter  that 
give  the  specified  output  for  the  given  input.  In  each  case,  the  problem  can 
be  converted  into  a  generalized  eigenproblem. 
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A  1000  X  1000  eigenproblem  requires  approximately  8  *  (1000)^  floating 
point  operations,  which  can  be  done  in  approximately  two  minutes  on  a  100 
MFlop  computer.  Several  sets  of  exponents  that  result  in  d  <  1000  solutions 
are  given  in  the  table  below.  The  number  of  floating  point  operations  required 
for  the  entire  solution  is  10[min(ni)d]®  which  is  given  in  the  last  column.  A 
100  MFlop  machine  does  10®  floating  point  operations  per  second,  so  the 
numbers  in  the  last  column  are  all  with  that  exponent. 


n  =  E?=i «» 

ni 

n2 

nz 

714 

ns 

ne 

j _  (ni+n2+...H-nib)J 

^  ni  !n2!...ni.l 

10[mm(ni)d]® 

6 

■1 

1 

1 

1 

a 

ni 

720 

37e-h08 

7 

2 

2 

1 

630 

256+08 

7 

3 

2 

2 

210 

7e+08 

7 

3 

3 

1 

140 

.3c+08 

8 

4 

4 

70 

2g+08 

32 

30 

1 

1 

992 

98e-f08 

1000 

999 

1 

1000 

1006+08 

1.2  Segre  Embeddings 

Multi-affine  equations  on  x  P"*  x  ...  x  P”*  can  be  rewritten  as  linear 
equations  on  p(»»i+i)(»»2+i)...(nfc+i)-i  using  the  Segre  embedding. 

For  i  =  1  to  k,  let 


Let  vq  =  1.  For  j  =  1  to  k,  let 


Vj-i 


g  p(nx+l)("2+l)—("i+l)~l  (10) 


Hns^j-1 

So  Vk  is  a  mapping  from  F"'  x  F"*  x  ...  x  F”*  to  p("i+i)(”2+i)-(«*+i)-i. 


1.3  A  Formula  for  the  Number  of  Solutions 


Theorem 

The  degree  of  the  Segre  embedding  of  F”i  x  F”*  x  ...  x  F”* 

in  ^  is: 

j_  (ni+Tt2  +  -  +  nji;)! 

ni!n2!...njk! 

So  intersecting  F”‘  x  F"*  x  ...  x  F”‘  with  a  codimension  ni  +n2  +  ...+nfc 
hyper-plane  (ni+n2  +  ...+nfc  hnear  equations)  in  p{«i+i)("2+i)  -K+i)-i  gives 
d  =  pouii-  solutions. 

Th^proof  below  uses  techniques  from  intersection  theory  that  can  be 
found  in  [Pulton, 84]. 

Proof 

Let  0  be  the  Segre  embedding: 

0  :  F"!  X  F"=*  X  ...  X  F"*  ->  (12) 

The  hyper-plane  class  u  in  -  pulls  back  to: 

ai  -f  02  +  ...  +  Ofc  e  F2(F”‘  X  F"^  X  ...  X  F”'=)  (13) 

where  is  the  2"''  Cohomology,  and  Cj  is  the  hyper-plane  class  in  F"j  . 
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Look  at  the  coefficients  in  the  expansion  of: 

(aj  +  02  +  ...  +  e  /f2(ni+n*+...+n*)(pni  pn*  ^  ^  pn»^ 

This  expression  can  be  expanded  by  repeatedly  applying  the  binomial  ex¬ 
pansion.  Many  of  the  resulting  coefficients  are  zero  because  for  each  i  from 
1  to  k  we  get: 

=  0  (15) 

We  only  get  one  nonzero  coefficient  in  the  expansion: 

(oi  -h  02  + ...  +  = 

(ni  +  n2  +  ...  +  n]b)(oi)(’“+”2+-+"*-^Ha2  +  03  -h ...  +  a*)  + ... 

d(o;*^aJ*...ajJ*)  + ... 

The  only  term  in  the  expansion  where  each  Oj  has  exponent  less  than 
rii  +  1  is  the  one  where  each  Oj  is  taken  nj  times. 

(ai  -f-  02  +  ...  +  afc)(»*i+"2+-+»**)  =0  +  d(aj‘a5*...a2*)  -f-  0  -I- ...  -1-  0  (16) 


The  coefficient  on  this  term  is: 


(ni  -fn2  -I- ...  -fTifc)! 

ni!ni!...nfc! 


(17) 


Therefore  d  is  the  degree  of  the  Segre  embedding,  and  there  are  d  solutions 
to  a  set  of  ni  +  Ti2  +  ...  +  nk  affine  equations  on  this  space. 


1.4  Dilation  Procedure  for  Forming  an  Eigen- 
problem 

In  this  section,  we  use  the  root  count  d  and  its  factorization  to  construct  an 
algorithm  for  computing  the  solutions. 
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We  will  work  with  ni  +712  + ...  H-n*  multi-affine  equations  on  x  F"*  x 
...  X  F”*.  For  i  =  1  to  k,  let 

■  1  ■ 

Xi,2 

Xi=  .  G  F”*  (18) 


L  J 

The  ni  -f  ri2  +  ...  +  njb  multi-affine  equations  in  the  Xij  variables  can  be 
written  as: 


Avk{Xu  X2, ...,  Xfc)  =  0  €  (19) 

where  the  matrix  A  has  ni-|-n2+...+nfc  rows  and  (ni-f-l)(n2+l)...(njfc-|-l) 
columns.  The  vector  VkiXi,  X2,  ...,Xk)  lies  in  the  Segre  embedding  of  F"i  x 
F”*  X  ..  X  F”*  in 

The  expression  Avk{XuX2,...,Xk)  =  0  has  fewer  equations  than  mono¬ 
mials.  We  need  to  append  more  equations  to  get  as  many  equations  as 
monomials.  We  can  multiply  the  equations  by  enough  monomials  to  get  as 
many  equations  as  monomials  in  the  new  matrix  equation.  Let 

sumXi  =  Id-  Xi^i  +  ...  -I-  (20) 

The  vector  Vfc(Xi,  X2, ...,  Xk)  contains  the  (ni-M)(n2+l)...(nfc+l)  monomials 
that  appear  in  the  expression: 

{sumXi){sumX2)...(sumXk)  (21) 


In  order  to  get  an  eigenproblem  with  d  solutions,  where 

\  nj  J  ni!ni!...nfc! 

we  can  use  an  equation  of  the  form: 

0  =  MiXk)w(Xi,X2,...Xk-i) 


(22) 

(23) 
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where  w(Xi,  X2,  ••^Xk-i)  contains  d  monomials.  To  get  these  d  monomials, 
we  proceed  as  follows. 

For  j  =  1  to  k,  let 

k 


Tj 

»=J 

(24) 

Then 

k-l 

<i=II 

j-i 

("t") 

(26) 

The  expression: 

(stX77xXj)’’^+‘ 

contains 

monomials 

(26) 

Let  the  vector  w(X)  e  contain  the  d  monomials  in  the  expression 

(sumXi)’’*  (51x771X2)*'®. ..(snmX/b-i)’’*  (27) 

To  get  these  d  monomials,  we  can  multiply  the  original  (tii  +  l)(n2  + 
l)...(77jfc  +  1)  monomials  contained  in  ujfc(X)  by  D  other  monomials  where 

£)  =  JJ  f  ^  —  rikd/ {n\  +  7x2  + ...  +  Tijk)  (28) 

i=i  \  / 

Let  g{X)  e  contain  the  D  monomials  in  the  following  expression: 

(sTXTnXi  {sumXi )  ^''®  .  (sumXk- 1 )  ^  (29) 


Let  G(Xi,  X2, ...,  Xk-i)  be  a  matrix  containing  each  of  the  D  monomials 
in  ^(Xi,X2,  ...,Xjb_i)  times  an  identity  matrix  with  tii  +  7x2  +  ...  +  tx*  rows. 


G'(Xi,X2,...,Xfe_i) 


*fni+n2+...  *nMX) 

*^ni  4-112+...  +nfcfl^(X) 


L  •^ni+n2+...  +nk9D{X)  J 


(30) 
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The  matrix  CiX^X^,  ...,Xk-i)  has  (ni  +  n2  +  ...  +  nk)D  =  n^d  rows  and 
(ni  +  n2  + ...  +  njk)  columns. 

The  following  then  are  rikd  equations,  in  {uk  +  l)d  monomials. 

. Xt.i)Avt(XuX2,...,Xu)  =  0  e  (31) 

The  above  expression  contains  the  {uk  +  l)d  monomials  found  in: 

{sumXiY^  {sumX2y^  ...{sumXk-iy‘‘  (sumXk)  (32) 

Separate  all  d  monomials  in  {Xi,X2,  ..■,Xk-i)  into  the  vector  w{X),  and 
all  n*  +  1  monomials  in  Xk  into  a  Ukd  x  d  matrix  M{Xk)- 

0  =  MiXk)w{Xu  X2, ...,  Xfc_i)  (33) 

The  Ukd  X  d  matrix  M{Xk)  is  affine  in  Xk 

M  (Xk)  =  Mk,o  +  Mk,iXk,i  +  ...  +  Mk,nk^k,nk  (34) 

The  Mk,i  matrices  are  linear  in  the  original  data  matrix  A. 

To  convert  this  system  into  an  eigenproblem  with  one  of  the  variables  as 
the  eigenvalue,  define  the  following  two  n^d  x  rikd  matrices: 

Mo  =  [Mfe,o,  Mk,i, ...,  Mk,nk-i]  (35) 

Mi  =  [Mk,nk,0,-,0]  (36) 

and  the  following  size  Ukd  vector: 

w 

Xk,lW 

Xk2fJJ 

w  =  ’.  (37) 


L  ^k,nk-l'^  J 

Then  the  following  equation  is  a  size  rikd  x  Ukd  generalized  eigenproblem 
with  at  most  d  finite  eigenvalues. 

0  =  M(Xk)w  =  (Mo  +  Xk,nkMi)  w  (38) 
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1.5  Deflation  Removes  Extraneous  Roots 

Generalized  eigenproblems  can  be  converted  to  regular  eigenproblems  if  ei¬ 
ther  Mo  or  Ml  is  full  rank.  Some  robotics  problems  have  roots  at  Xk,nk  = 
which  are  extraneous.  The  roots  at  ~  0>  removed  by 

deflation. 

If  Mq  is  low  rank,  then  there  are  roots  at  Xk,nk  —  0- 

If  Ml  is  low  rank,  then  there  are  roots  at  Xk,nk  =  oo. 

If  Mq  +  iMi  is  low  rank,  then  there  are  roots  at  Xk,nk  ~  ■*"*• 

If  Mo  —  iMi  is  low  rank,  then  there  are  roots  at  Xk,nk  =  ~*- 

If  Mo  +  xMi  is  low  rank  for  generic  x,  then  the  solution  set  contains  sets 
of  positive  dimension,  rather  than  just  points. 

The  rest  of  this  section  uses  deflation  to  eliminate  roots  at  infinity.  By 
interchanging  the  roles  of  Mo  and  Mi,  the  roots  at  0  can  also  be  eliminated. 
To  remove  roots  at  ±i,  first  do  a  linear  fractional  transformation  that  maps 
-fi  to  0  and  maps  —i  to  oo. 

If  the  system  has  repeated  Jordan  blocks  with  roots  at  the  point  being 
eliminated,  the  deflation  procedure  reduces  the  size  of  each  of  those  Jordan 
blocks  by  one,  each  pass  through  the  loop.  So  the  maximum  number  of 
passes  is  equal  to  the  maximum  size  of  any  Jordan  block  associated  with  the 
eigenvalue  being  eliminated. 

The  matrices  Mo  and  Mi  in  the  last  section  are  square,  rikd  x  Ukd,  but 
the  Ml  matrix  multiplying  Xk,nk  Is  low  rank.  (For  matrices  with  more  rows 
than  columns,  see  the  section  on  "Incomplete  Intersection”.)  To  get  rid  of 
the  roots  at  infinity,  we  can  get  a  smaller  d  x  d  generalized  eigenproblem 
using  deflation. 
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THEOREM  (Eigen-Problem  Deflation) 

The  finite  pair  A,  e  is  a  solution  to  0  =  \Mq  +  AJlfiJe  iff  the  finite 
pair  A,  /  is  a  solution  to  0  =  [Lo  +  where  Lq  and  Li  are  defined 

in  terms  of  the  svd  of  the  Mi  as  follows: 


Mx  =  [UuU2] 


El  0 
0  0 


y* 


(39) 


Then  taking  the  svd  of  U2*Mq‘. 


U2*Mq  =  [UiM 


Si  0 
0  0 


[Vi.Va]* 


(40) 


Using  these  svd’s,  Lq  and  Li  are  defined  by: 

[Lo  +  ALi]  =  Ui*[Mo  +  AMilVa  (41) 

A  solution  pair  A,  /  of  0  =  [Lq  +  ALi]/  is  then  mapped  back  to  a 
solution  pair  A,  e  of  0  =  [Mo  +  AMi]e  by: 


e  =  V2/ 


(42) 


PROOF 

case  1)  =  0  and  V2  7^  0  (infinite  number  of  finite  A,  e  solutions) 

This  can  only  happen  if  Mi  =  0.  In  this  case,  Ui  is  the  empty  set,  so  the 
equation  [Lq  +  ALi]/  =  0  is  satisfied  by  any  A  and  any  /.  Since  Mi  =  0,  U2 
is  square,  so  V2  =  kernel{Mo)  =  kernel{Mo  +  AMi),  so  e  =  V2/  is  a  solution 
to  [Mo  +  AMije  for  any  A  and  any  /. 

case  2)  V2  =  0  (no  finite  A,  e  solutions) 

This  can  only  happen  if  U2*Mo  is  full  column  rank.  In  this  case,  1/2* [Mo  + 
AMi]  =  U2*Mo  which  has  no  right  kernel,  so  [Mb  +  AMi]e  cannot  be  zero 
for  any  nonzero  value  of  e.  In  this  case,  there  are  no  nonzero  solutions  for  e. 
Also,  since  V2  =  0,  /  =  0-  So  in  this  case,  neither  system  has  a  solution. 

case  3)  Ui  ^  0  and  V2  ^  0  (finite  number  of  finite  A,  e  solutions) 
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Let 


(43) 


{  =lVi,V2re 
*  [  J  . 

Multiplying  the  eigenproblem  by  the  square  orthogonal  matrix 
gives: 


0  =  [Ui,  t/a]* [Mo  +  AMi]e  =  [Ui,  C/a]* [Mo  +  AMi][Vi,  Va] 

(44) 

Let 

1,1,1  =  t/nMo  +  AMi]Vi 

(45) 

Li,a  =  U;[Mo  +  AMijVa 

(46) 

1^2.1  =  t^2*[Mo  +  AMi]Vi 

(47) 

1^2.2  =  C4lMo  +  AMi]Va 

(48) 

Then:  _ 

(£/„I/J'|J4o  +  AMJ[Vi,V2]- 

(49) 

Since  U^Mi  =  0  and  C/jMoVa  =  0,  the  La, 2  block  is  identically  zero. 

With  La, a  identically  zero,  the  bottom  part  of  the  matrix  equation  gives: 

0  =  La.i/  +  La, a/  —  L2,i/ 

(50) 

Plugging  in  the  expression  for  La,i  then  gives 

0  =  U;[Mo  +  AMi]Vi/  =  iU*Mo)Vj  =  UiSj 

(51) 

Since  the  matrix  UiSi  has  full  column  rank,  this  implies  that 

o 

II 

(52) 

o 

II 

(53) 
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This  reduces  the  eigenproblem  to: 

0  =  Li.2/  (54) 

which  can  also  be  written  as: 

0=[C/nMo  +  AMi]V2]/  =  [Lo  +  ALi]/  (55) 


After  solving  this  new  problem,  the  vector  e  can  then  be  obtained  from: 

e  =  V2/  (56) 

The  above  procedure  will  have  to  be  repeated  if  the  resulting  L\  matrix 
is  still  low  rank  (replace  Mo,  Mi,  and  e  with  Lo,  Li,  and  /  then  repeat).  The 
eigenvectors  get  mapped  back  by  the  product  of  the  V2  matrices  from  each 
step  of  the  iteration. 

After  the  final  Li  is  full  rank,  the  following  dxd  eigenproblem  gives  the 
d  values  of  A  =  xjfc.n*  and  the  corresponding  d  values  of  the  eigenvector  / . 

[Lo  +  Xk,nk^i]f  =  0  (57) 

The  vector  w{Xi,  X2, Xk-i)  is  the  first  d  elements  of  the  vector  w.  The 
other  Xij  can  be  extracted  from  the  vector  w{Xi,X2,  ...,Xjfc_i). 

1.6  Pseudo  Code  for  the  Dilation  Procedure 

Let 

M=[Mo,Mi,...,M„J  (58) 

Let  _  _ 

wiX) 

Xk,iw(X) 

Xk,2‘w{X) 

w{X)  =  .  (59) 

Xk;nkW(X) 
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Then 


Mw{X)  =  G{X)Avk{X)  (60) 

When  multiplying  monomials  in  G{X)Avk{X)  to  get  Mw{X),  the  expo¬ 
nents  on  the  monomials  in  G(X)  and  njk(X)  add.  Let  basea{-,j)  be  a  column 
of  exponents  on  the  Xi,m  in  the  monomial  in  G(X).  Let  the  function  index 

be  defined  such  that  indexG{basea(}ij))  =  j- 
The  monomial  in  C?(X)  is  given  by: 

n  n  (61) 

J=1  m=l 

The  r*'*  monomial  in  Uifc(X)  is  given  by: 

n  li  (62) 

i=l  m=l 

The  product  of  the  monomial  in  G{X)  and  the  monomial  in  ^^(X) 
gives  monomial  number  s  =  index^ibasegi:,])  -I-  basev(:,Tn))  in  u)(X).  The 
monomial  in  w(X)  is  given  by: 


TT  TT 

11  11 

1=1  m=l 


(63) 


These  base  and  index  ideas  can  be  used  in  the  following  pseudo  code  for 
filling  in  the  M  matrix. 

for  i  =  1  to  D  Put  iV  =  ni  -1-  712  -I- ...  +  rows  of  A  into  M 

for  j  =  1  to  (ni  +  l)(7i2  +  l)..-(wfc  -f  1)  Put  a  column  of  A  into  M 

Compute  exponents  baseG{:,i)  on  the  i^  monomial  in  G{X) 
Compute  exponents  base„{:,j)  on  the  monomial  in  Vk(X) 
row^  =  index  ^{baseG{';i)  +  base„{:,j)) 
rowsfy  =  (i  —  l)*N+l:i*N 
M(rowSj^,row^)  =  A{:,j) 
end 

end 
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A:  -  1  copies  of  1  copy  of  P”*. 


1.7  Examples 

P^  X  ...  X  P^  X  P"*  X  P^ 

Get  dy.  d  matrix  Mo  +  XkMi. 

{sumXiY^  (sumX^y^  ...{sumXk-iY'‘  (sumXk) 

=  (1  +  a;i)"*+*-=*(l  +  +  a;*-!,!  +  ...  +  Xk-i,m)\l  +  Xk) 

So  M{xk)  has  two  monomials,  1  and  Xk,  while 

w  has  d  =  (m  +  (A:  —  l))(Tn  +  (A:  —  2))...(m  +1)  =  (m  +  (A:  —  l))!/m! 
monomials. 


In  particular,  P^  x  P*”  x  P^ 

(sumXiY^  {sumXk-iY^  (sumXs) 

=  (1  +  +  2:2,1  +  ...  +  a:2.m)Hl  +  2:3) 

Let[l,Xl]eP^ 

r  1 

2:2,1 

2:2,2 


gpm 


^2,m 


and  [1,2:3]  e  P^. 

The  2  ♦  (m  +  1)  *  2  monomials  form  the  vector  w: 


m  = 


il 

Xi3L 

xz3i 


X1X3U  J 


(64) 


(65) 


The  2  +  m  equations  in  2  ♦  (m  +  1)  ♦  2  monomials  are  in  the  form: 


0  =  [v4oo,.<4oi,  Aio,.i4ii]u; 
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Dilation  is  done  by  multiplying  the  above  equations  by  all  monomials  in 
(1  +  a:i)’"+S  giving: 


v 

XiU 

X\^V 


■Aqo 

-<4oi 

0 

0 

A\q 

An 

0 

. 

0  ■ 

0 

-^00 

•^01 

0 

• 

• 

0 

-^10 

0 

• 

0 

• 

0 

•^00 

Aqi 

0 

• 

• 

0 

-^10 

^11 

0 

0 

, 

0 

Aoa 

.^1 

0 

. 

• 

0 

-^10 

^11 . 

XzV 

X3X1V 

X3Xi^2L 


This  equation  can  be  written  in  the  form 


(67) 


0  =  [Mo  +  xaMjuec 


(68) 


where  Md  and  Mi  are  each  {m  +  2)  (m  +  1)  x  (m  +  2)  (rn  +  1) .  This  gives  a 
square  eigenproblem  with  (m  +  2)(m  +  1)  eigenvalues. 
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Chapter  2 

More  Equations  than 
Unknowns 


2.1  Curves  on  Segre  Spaces 

Let  a  be  the  Segre  map  from  (P^)"  x  P"*  to  . 

a  :  (P^)"  xP”^  (1) 

Let  Y  C  P^  be  a  linear  subspace  of  P^. 

Let  C  be  a  degree  d  curve  given  by 

c  =  rnE  (2) 

where  the  intersection  would  be  empty  if  Y  were  in  general  position,  so  the 
intersection  of  Y  and  E  is  not  transverse. 

Let  a'  be  the  Segre  map  from  (P^)"^  x  P"*  to  P^' .  Let  Y'  C  P^  be  a 
linear  subspace  of  P^' . 

QUESTION:  What  is  the  smallest  n'  such  that 

C  U  {points}  =  1^*0  2'  (3) 

where  some  extraneouse  points  have  been  introduced  in  the  elimination  of 
some  variables. 

ANSWER(?):  n'  =  smallest  I  s.t.  d<d!  =  {I  +  m)\/m\ 
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SOLVE:  Use  dilation  to  form  &df  xd!  parameterized  eigenproblem.  Use 
deflation  to  eliminate  the  extraneous  points,  leaving  a.  d  x  d  parameterized 
eigenproblem. 


2.2  Example:  7R  mechanism 


For  each  fixed  value  of  zj,  we  get  30  +  m  equations  on  x  x  x 
P^  X  P^  X  P^,  with  16  solutions.  If  only  6  equations  are  kept,  get  6!  =  720 
solutions. 

If  30  or  more  equations  are  kept,  a  linear  algebra  procedure  can  be  used 
to  eliminate  some  variables  and  get  a  smaller  problem  with  only  the  16  le¬ 
gitimate  solutions.  However,  if  this  elimination  procedure  is  taken  too  far, 
it  ends  up  as  4  equations  on  P^  x  P^  x  P^  x  P^  with  4!  =  24  solutions. 
Eight  of  these  24  solutions  are  at  ±i  which  are  extraneous.  The  remaining 
16  solutions  are  legitimate. 

To  see  how  the  extra  equations  eliminate  the  extraneous  solutions,  we 
can  look  at  the  Smith  form  of  the  affine  matrices.  The  30  -I-  m  equations  are 
of  the  form: 


0  =  [Mrfl  -h  27Af7,i] 


h2 

■fl6 

h 

Z&Iz2 

zJz  . 

V%{Zz,Z2,Zi) 


For  each  fixed  value  of  zt,  let 


=  [Afe.o  +  =  [-^7,0  + 


I32 


(4) 

(5) 


The  matrix  Mq{zq)  is  size  (30+m)  x  32.  It  is  rank  30  for  all  zq  (unless  m  = 
0,  then  rank  drops  for  discrete  values  of  zg)-  Its  Smith  form  decomposition 
is: 


^5(26)  = 


^6top 


where  the  size  (30  +  m)  x  30  matrix  Geiz^)  =  [G'6£,(«6),C?6|i]  is  rank  30 
for  all  zg  (unless  m=0,  then  rank  drops  for  discrete  values  of  zg).  The  matrix 
iiCg(zg)  is  size  30  x  32.  The  top  16  rows  of  K^ize)  are  constant,  denoted  by 
the  16  X  32  matrix 
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(7) 


Since  Geize)  is  full  rank  for  all  ze, 

kernel[Me{z6)]  =  kernel[K6(ze)] 

In  this  case  we  also  get  (WHY?): 

kernellKeize)]  =  kernel[K6top]  (S) 


so  we  get  a  new  problem,  with  ze  eliminated,  and  only  16  equations: 


0  =  ^etop 


he' 

Zhhe  _ 

Zih  _ 

V8iZz,Z2,Zi) 


(9) 


If  we  repeat  the  above  procedure  again,  eliminating  2:5,  we  introduce  some 
extraneous  roots  because  this  time  the  kernel  of  the  new  K{z)  matrix  is  not 
the  same  as  the  kernel  of  the  constant  part  of  K{z). 


Let 


Msizs)  =  K^top 


he 

zzhe 


(10) 


The  16  X  16  matrix  M^iz^)  is  rank  12  for  all  (generic?)  yalues  of  zs,  and 
can  be  factored  as: 


Msizs)  =  lG5x,(^5)><J5fl] 


Kziop 


(11) 


where  the  size  16  x  12  matrix  ^5(25)  =  [05^(25),  ^5^]  is  rank  12  for  all 
(generic?)  values  of  zs.  The  matrix  ^£^5(25)  is  size  12  x  16.  The  top  4  rows 
of  K^iz^)  are  constant,  denoted  by  the  4  x  16  matrix  K^top- 


The  commutative  diagram  below  shows  what  happens  if  we  keep  only  the 
constant  part  oi  K{z)  at  each  step.  The  polynomial  kernels,  X{z),  of  M{z) 
are  represented  as  coefficient  matrices,  X(bj)  —  times  polynomials 

Pj{z)  using  Lagrange  interpolation.  Lagrange  interpolation  results  in  X(bj) 
being  the  same  rank  for  each  j,  while  the  standard  polynomial  basis,  X(z)  = 
Xiz\  gives  coefficient  matrices  that  can  have  lower  rank  for  Xq  and/or 
Xi^^.  Lagrange  interpolation  also  allows  us  to  compute  only  numerical 
kernels  of  M{bj)  rather  than  using  Wedderburn  theory  to  compute  X(z)  = 
M{z)^. 
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After  the  following  diagram  stabilizes,  the  final  M{z)  has  as  many  equa¬ 
tions  as  unknowns.  All  original  roots  survive,  but  some  ’’extraneous”  roots 
may  be  obtained.  When  the  final  M{z)  has  more  columns  than  rows,  we  can 
use  dilation  to  produce  a  square  eigenproblem.  The  "extranous”  roots  are 
either  at  ±i  when  z  =  tan{6/2)  or  at  0,inf  when  z  =  eU0).  The  possible 
reason  for  the  ’’extraneous”  roots  is: 

The  discarded  non-constant  part  of  K(z)  would  have  removed  those  roots. 
Recall  that  the  real  part  of  the  quaternion  equations  were  dropped,  because 
they  were  not  multi-afiine.  Since  we  know  where  the  extraneouse  roots  are, 
we  can  remove  them  using  deflation,  before  an  eigenproblem  is  done.  When 
Wedderbum  gave  X(z)  with  low  degree,  it  was  in  eflfect  removing  roots  at 
infinity,  similar  to  the  way  deflation  works. 

The  second  commutative  diagram  uses  Wedderbum  theory  to  explicitly 
constmct  X(z)  =  M^z)-^.  The  degree  of  the  X{z)  polynomial  can  be  less 
than  its  generic  value  if  the  equations  0  =  M{z)v  has  ’’extraneous”  roots 
at  infinity  or  zero.  Roots  at  infinity  or  zero  can  cause  Xi  =  0  for  the  first 
few  or  last  few  values  of  i.  In  either  case,  replaceing  X  (z)  with  the  lower 
degree  pol5momial  obtained  by  dropping  the  zero  valued  coefficients,  gives  a 
smaller  system  of  equations  which  no  longer  contain  the  ’’extraneous”  roots 
at  infinity  or  zero. 

The  advantage  of  using  the  Lagrange  interpolation  instead  of  the  Wed¬ 
derbum  theory,  is  that  the  degree  of  the  X{z)  can  be  "fractonal”,  ie  the 
leading  coefficient  matrix  could  be  low  rank.  This  causes  problems  when 
computing  Mn_i(zn-i)  using  Wedderbum.  However,  Lagrange  interpolation 
has  no  such  problems. 

When  the  final  (or  original)  M{z)  is  square  (eg  3P-3R  stmcture),  then 
only  some  roots  are  legitimate,  while  many  roots  are  ’’extranous”.  Some  of 
these  are  the  ±i  or  0,  inf  roots  known  to  satisfy  the  multi-affine  equations, 
while  others  do  not  have  eigenvectors  that  satisfy  the  required  Segre  stmc¬ 
ture.  For  the  3P-3R  structure,  this  can  be  avoided  by  using  a  smaller  system 
of  equations  associated  with  only  the  rotation  equations  (no  translation  equa¬ 
tions). 

First  commutative  diagram  for  the  6R  structure. 
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0 


0 


0 


0 


Pi{z^)h 

PoizA)h 


R}^ 

P2{z5)h 

Pi(25)/4 

L  Poiz5)h 


R^^ 

P7{.Ze)h 

L  PQi.Z^)h 
R^ 


[X,ih),X,{bo)] 


X,(z,) 


R^  0 

th 

M,{z,) 

R  — >  ^  0 

h 


[Xe(b7),...,Xe(bo)] 


Z4I8 

R}^ 


the 


he 

Z5I16 


[M40,  M41] 


i/4 


0 


M^iz^) 


R}^  ->  R'^ 

■1'  ■^16 


32 


[M50,  M51] 


X6(Z6) 


R 


t/32 


R}^  0 


MeCze) 


paa+m  0 

(12) 


Alternatively,  the  8  x  8  quadratic  matrix  M^i^z^)  =  [-X8^^{z5),X8^X^8)\ 
gives  an  eigenproblem  with  16  roots. 

Second  commutative  diagram  for  the  6R  structure. 
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zJi 

h 


^5top  (■^5) 
•^Skot  ('^5) 


M4{z5) 


M^izs) 

Rio  RlO 

4-  [—^ihth] 


0 


(13) 


For  other  problems,  such  as  the  4R-2P  structure,  we  can  first  compute 
■^4(24),  then  apply  the  above  procedure  twice  to  obtain  a  2  x  2  degree  4 
matrix  with  8  good  eigenvalues.  Generic  data  would  result  in  a  2  x  2  degree 
6  matrix  with  12  eigenvalues.  If  we  kept  the  zero  coefficients  in  the  2x2 
degree  6  matrix,  we  get  8  "good”  eigenvalues,  and  4  at  zero  or  infinity. 


2.3  Incomplete  Intersections  (for  eigenprob- 
lems) 

After  a  generic  linear  change  of  coordinates,  any  0-dimensional  ideal  in  n  vari¬ 
ables,  containing  d  non-repeated  points,  can  be  represented  by  n  polynomial 
equations  (a  complete  intersection).  The  roots  of  this  set  of  n  polynomi¬ 
als  will  be  the  original  d  points  (plus  additional  points  at  infinity  in  the 
homogeneous  case).  However,  if  some  structure  is  imposed  on  the  type  of 
polynomial  equations  used  to  represent  the  ideal,  e.g.  multi-afiine  pol3mo- 
mials,  then  it  niay  be  necessary  to  use  more  equations  than  unknowns  (an 
incomplete  intersection). 

Given  any  system  of  polynomial  equations  of  any  degree,  it  can  be  con¬ 
verted  to  a  system  of  multi-affine  equations  by  replacing  powers  of  variables 
with  new  variables.  This  generates  new  variables  as  well  as  new  relation¬ 
ships  amongst  the  variables.  If  the  ideal  is  O-dimensional  with  d  non-repeated 
point  solutions,  then  after  a  generic  linear  change  of  coordinates,  the  reduced 
Grobner  basis,  G  =  {gi,9u of  the  ideal  is  of  the  form  [Gianni  and 
Mora  1989]: 


21 


Xi 

r  ^0  1 

’  9i  ’ 

X2 

92 

• 

= 

-B 

• 

• 

• 

• 

Xn-l 

-  9n  . 

X^ 

L  J 

J 

where  the  matrix  B  is  unique. 


(14) 


Since  d  can  be  any  integer,  it  need  not  be  of  the  form  obtained  for  generic 
multi-affine  systems,  in  particular  it  need  not  even  be  factorable.  When  a 
system  of  polynomials  of  any  degree  is  put  into  multi-affine  form,  coefficients 
on  some  monomials  may  be  missing.  This  may  result  in  some  roots  being  at 
infinity  or  more  equations  than  unknowns  (incomplete  intersection).  Some 
systems  of  polynomials  for  0-dimensional  ideals  may  also  start  out  as  incom¬ 
plete  intersections.  When  there  are  more  equations  than  unknowns,  then  the 
procedure  described  in  previous  sections  multiplies  them  all  by  new  mono¬ 
mials,  resulting  in  an  ’’eigenproblem”  with  more  equations  than  unknowns. 
These  extra  equations  eliminate  some  of  the  eigenvalues  and  eigenvectors. 


An  m  X  n  eigenproblem  with  d  solutions  has  d  eigenvalues  and  d  eigen¬ 
vectors  (assume  all  Jordon  blocks  are  size  1?).  Let  C  and  D  be  the  d  x  d 
diagonal  matrices  whose  ratio  give  the  d  generalized  eigenvalues,  and  let  Vo 
be  the  nxd  matrix  whose  column  is  the  eigenvector.  The  eigen  problem 
can  be  written  as; 

0  =  [Ao,  Bo] 

where  Aq  and  Bo  are  each  mxn  matrices  with  m>n>  d. 


VoC 

-VoD 


(15) 


We  want  to  compute  Vo,  C,  and  D. 

If  we  numerically  compute  the  kernel  of  the  [Ao,  Bo]  matrix,  we  get 


Bi 

-Ai 


kernel[Ao,  Bo 


(16) 


then 


VoC 

Bx' 

-VoD 

.  -^1  . 

(17) 
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where  Vi  is  an  unknown  di  x  di  matrix, 

Since  the  diagonal  C  and  D  matrices  commute,  right  multiplication  of  the 
top  half  of  the  above  equation  by  D  gives  the  negative  of  right  multiplication 
of  the  bottom  half  of  the  above  equation  by  C.  This  gives  a  new  n  x  di 
eigenproblem: 

0  =  IA„S,1  (18) 

Repeating  this  process  gives: 

=  kernel[Ai,  Bi]  (19) 

—A2 


then 


II 

0  Q 

1 

■^2  1  y 

(20) 

where  is  an  unknown  ^2  x  da  matrix, 

iProm  this,  we  get  a  new  di  x  d2  eigenproblem: 

0  =  [A2,  B2] 

V2C  ■ 

-V2D 

(21) 

This  procedure  is  repeated  until  the  integer  sequence  dj  converges  to  d. 
If  the  [>lo,  Bo]  matrix  starts  out  square,  then  numerical  experiments  indicate 
that  the  integer  sequence  converges  for  i  >  ij  where  ij  is  the  size  of  the 
largest  Jordan  block  (when  0  <  d  <  00  ). 

The  eigenvectors  of  the  original  system  can  be  recovered  using  either: 

VbD‘  =  AiViD^-^  = 

...  =  AiA2...AiVi 

(22) 

or 

VoC^  =  BiVi(7-^  = 

...  =  BiB2—BiVi 

(23) 

The  Hiagnnfll  entries  of  the  C  and  D  matrices  are  never  both  zero  at  the 
same  location  on  the  diagonal.  Therefore  any  column  of  the  Vq  matrix  can 
be  solved  for  using  one  or  the  other  of  the  above  two  equations.  If  Djj  ^  0, 
then 

colj[Vo]  =  [A,A.,...A,]cdi[VAIDjf  (24) 
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If  ^  0,  then 


colj[Vo]  =  [BiB2...Bi]colj[Vi\/Cj/  (25) 

Actually,  the  division  by  the  diagonal  entries  oi  C  or  D  can  be  eliminated, 
since  eigenvectors  are  typically  scaled  to  unit  length  anyway. 

Note  that  if  we  multiply  the  equation: 

VqD^  =  AiA2...AiVi  (26) 

on  the  right  by  we  get  the  same  thing  as  multiplying  the  equation 

Vo&  =  BiB2...BiVi  (27) 

on  the  right  by  D'.  Therefore  we  get: 

r  V-O*  1 

0  =  [AiA2...A<,  BiB2...Bi]  (28) 


CONVERGENCE  OF  THE  ITERATION 


THEOREM:  If  the  m  x  n  eigenproblem  [Ao  +  XBo]v  =  0has0<d<oo 
solutions,  then  the  iteration  defined  by: 


0 = [A,Bi]  f 


converges  in  imax  +  1  iterations  to  a  square  d  x  d  eigenproblem,  where  imax 
is  bounded  by: 

ceiling  [— — — ]  <  imax  <n  —  d  (30) 


m  —  n 


See  Gantmacher  Normal  Form  to  get  an  exact  equation  for  imax  [Thomps,1970]. 
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PROOF 

Let  d-i  =  m,  do  =  n,  and  for  i  >  0,  let  dj+i  be  the  dimension  of  the 
kernel  of  Mi.  The  matrix  Mi  has  2di  columns,  so 

di+i  =  2di  —  rank{Mi)  >  0  (31) 


Note  that 


rank{[Ao  +  ABo])  =  rank  ^[i4o,  Bq]  ^  ^  <  rank{[Ao,  Bo])  (32) 


Assume  rank(^Mo)  ^  do  so  we  get  at  most  a  finite  number  of  solutions  to 
the  original  eigenproblem.  We  must  have  rank{Mi)  >  di,  else  we  would  get 
an  infinite  number  of  solutions  to  the  eigenproblem,  0  =  {Ai  +  ABi)u(A) 
for  all  values  of  A.  This  would  map  down  to  an  infinite  number  of  solutions 
to  the  previous  eigenproblem,  0  =  (Aj_i  +  ABi_i)Biu(A)  for  all  values  of  A. 
Eventually,  this  would  map  down  to  an  infinite  number  of  solutions  to  the 
original  eigenproblem,  0  =  (Ao  +  ABo)BiB2...Biu(A)  for  all  values  of  A. 

Because  rank{Mi)  >  di,  we  get: 

di+i  =  2di  —  rank{Mi)  <  di  (33) 


so  the  sequence  of  nonnegative  integers,  dj,  is  monotonically  decreasing. 

LEMMA: 

The  integer  sequence  dj  converges  to  d,  the  number  of  solutions  to  the 
original  eigenproblem. 

PROOF  OF  LEMMA: 

If  di^^+i  =  di^axi  we  get  a  square  eigenproblem  at  that  step.  If  the 
eigenvalues  are  not  repeated,  then  there  is  a  full  set  of  di^^,  eigenvectors, 
so  the  kernel  of  the  2dj^^2  x  di^^^  matrix  [Ai„„+i,Bi„„,+i]  will  be  di^„ 
dimensional,  so  dj  =  di„„^  for  all  i  >  imax- 

end  of  PROOF  OF  LEMMA: 

If  di  decreases  at  the  slowest  possible  rate,  i.e.  di+i  =  dj  — 1  until  i  —  imax^ 
then  imax  —n-d.  In  general 

imax  <n-d  (34) 
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Since  rank{Mi)  <  NumRows^Mi)  —  di-i,  the  equation 

di+i  =  2di  —  rank{Mi) 


(35) 


gives: 


rff+i  ^  di—\ 


Combining  this  with  the  initial  conditions,  d-\  =m,  do  =  n,  gives: 


(36) 


di>n  —  i{m  —  n) 


(37) 


This  decreasing  integer  sequence  converges  for  i  >  imax  where  imax  is 
obtained  by  setting  di  =  din  the  above  formula. 


imax  >  ceiling 


dfj  —  d 


=  ceiling 


m  —  doj 

Combining  the  upper  and  lower  bounds  gives: 

n  —  d 


n  —  d 


m  —  n 


ceiling 


m  —  n 


^  imax  —  ^  d 


(38) 


(39) 


To  obtain  a  square  eigenproblem,  one  more  iteration  is  needed,  since 
Ai  +  XBi  is  size  di_i  x  df. 


.dxd 


(40) 


When  the  sequence  converges  to  some  integer  d,  that  must  be  the  number 
of  solutions  to  the  original  eigenproblem  because  we  get  square  d  x  d  matrices 
Ai  and  Bi,  so  the  following  eigenproblem  0  =  [^4*  —  ABjju  has  d  eigenvalues. 
The  eigenvector  matrix  gets  mapped  back  into  the  original  large  space  by 
recursively  appl)dng  Cvi-i  =  BiVi  or  Dvi-i  =  AiVi. 

Ok  even  if  there  are  repeated  eigenvalues,  and  not  a  full  set  of  d  eigen¬ 
vectors. 
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COMMUTING  DIAGRAM: 


Let  Mi  =  [AjjBi].  Let  Li+i  =  kernel{Mi)  so  that  MiLi+i  =  0.  Define 
— Ai+i  and  Bj+i  as  the  bottom  and  top  halves  of  the  Lj+i  matrix. 


Bi+i 
— Ai+i 


=  L, 


'»+i 


so  AiBi+i  =  BiAi+i. 
Let 


0 

Ai 


and 


Bi  0 
0  Bi 


(41) 


(42) 


The  row  in  the  following  diagram  gives  the  eigenproblem.  The 
vertical  maps  show  how  the  eigenvectors,  kernels,  and  ranges  are  mapped 
(assuming  that  C  =  I).  The  following  diagram  commutes  since  AiBi+i  = 

BiAi-^i . 


J^dimax+2 


R^dir, 


C+1 


0 


0  — >■  b:^ 

B4 

0  — )■ 

IBs 

0  — y 

IBs 

0  — 


Ms 

— y 

iBs 

M2 

R2d^ 

— y 

iB2 

Ml 

R2d^ 

— y 

Mo 

R2dc 

— y 

R^^  — ^  0 

IB2 

R**!  — >■  0 

IB^ 

— >•  0 

i  Bo 

/r  — )•  0 


Note  that  if  the  assumption  C  =  /  is  replaced  with  the  assumption  D  =  /, 
then  a  similar  diagram  can  be  made  by  replacing  the  vertical  Bi  and  Bi  maps 
with  vertical  Ai  and  Ai  maps. 


If  di  =  d  for  i  >  ij,  then  we  get  the  following  d  x  d  eigenproblem  with 
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eigenvalues  =  diag(D).. 


0  — 


V 

-VD 


(44) 


After  the  previous  diagram  converges,  we  get  a  pull-back  from  the  follow¬ 
ing  commuting  diagram: 


Jpdij+i 

[Aij+i,Bij+i] 

•i'  ^ij+l 

-I'  ^ij 

(45) 

R^j 

j  -^ij] 

Rf^j-^ 

Pull  Back  =  ker  -[Aij,Bij]) :  R'^-'  0  ^  R"-'  (46) 


[4. 

0 

[4, 


0  ■ 
Bij+i 
~Aij+i  . 


(47) 


dim{PullBack)  =  (d.^  +  2dij)  -  rank{[Bij,  -(Ai^,  B^)])  (48) 

=  {dij  +  2dij)  -  rank{[Aij,Bij])  =  {dij  +  2dij)  -  {2dij  -  djj+i)  =  2d  (49) 
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2.4  Incomplete  Intersections  (general  poly¬ 
nomials) 


Any  set 
form  as: 


of  real  polynomial  equations  in  a:  €  C”  can  be  written  in  matrix 


0  = 


Mo 

Qoix) 


Vo{x) 


(50) 


Where  Mo  6  is  a  matrix  of  coefficients,  uo(a;)  6  Z^lx]  is  a  vector  of 

monomials,  and  Qo{x)  e  Z^'^[x]  represents  the  relations  among  the  monomi¬ 
als  voix). 


Qoi^)  —  Qo.t'^i(^) 


(51) 


i=l 


where  the  Qo.t  matrices  are  filled  with  integers,  and  the  Wi{x)  are  monomials. 
For  example,  in  the  generalized  eigenproblem,  0  =  [A  —  XB]eigvec,  we  get 

—eigvec 


Mo  =  [A,  B],  Voix)  = 


Xeigvec 


,  and  Qo(ic)  = 


In  the  general  case,  vo(x)  can  be  expressed  as; 

uo(x)  =  Mo-^co(x) 

for  any  polynomial  vector  Co(x)  that  satisfies 

Qo(x)Mo-‘-co(x)  =  0 

This  can  be  expanded  as: 

Jc 

0  =  Qo(:i^)Mo-^co(x)  =  ^Qo,m(x)Mo-^co(x) 

1=1 


=  QOy^Mo^,  — ,  Qo.kMo'^] 


Wiix)coix) 

W2ix)coix) 


Let 


[  u;ik(a;)co(a:)  J 

Ml  =  [Qo,iMo^,Qo,2MQ'^,...,QojiMo^] 


(52) 


(53) 


(54) 


(55) 


(56) 
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r  wi(a;)cb(x)  ■ 

,  .  W2(x)co(x) 

vi[x)  =  . 

_  Wk{x)coix)  . 

We  then  have  a  new  matrix  equation 

0  =  MiVi{x) 

The  new  monomial  vector,  Vi{x)  satisfies  the  relations: 

0  =  Qi(a;)ui(a;) 


where 


Qi{x)  = 


W2{x)I  —Wi{x)I  0  0 

wslx)!  0  -wi(x)I  0 


Wk{x)I  0 


0  0  ...  —wi{x)I 


This  gives  the  new  system: 


Ml 

Qiix) 


vi{x) 


If  we  can  solve  this  new  system  for  then  we  can  reconstruct  uo(2;) 

as  follows: 

uo(a;)  =  Ma^Cdix)  =  MQ^PiVi{x)lwi{x)  (62) 

where 

Pi  =  [0,0,...,0,/,0,...,0]  (63) 

such  that  PiVQ{x)  =  Wi{x)co{x). 

If  we  cannot  solve  the  new  system  for  Ui(x),  then  repeat  the  above  process 
until  Qj{x)Mj  converges  to  a  fixed  dimension  and  degree. 

Given  the  set  of  d  solutions:  X  =  [x(l),a:(2),  ...,a;(d)] 

Let  Vo  =  K(x(l)),vo(a;(2)),...,vo(x(d))].  Then  the  above  matrix  equa¬ 
tion  can  be  written  as: 
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0  =  MqVq 

(64) 

We  can  express  Vo  in  terms  of  the  numerical  kernel.  Mo'*', 

of  the  Mo 

matrix. 

VJ,  =  Mo-^Go 

(65) 

where 

Co  =  [co(x(l)),Co(a;(2)), ...,  Co(a;(d))] 

(66) 

Special  Case 
For  the  eigen-problem, 

1 

Xi 
XT. 

Xn-l 

This  gives  Qq{x)  =  [x„I„, 

In  a  slightly  more  general  case 


Vo  (a;)  = 


In 

Xjiln 


Vo  (a;)  = 


In 

Gq{x) 


1 

Xi 

XT 


Xn-1 


(67) 


(68) 


where  Go(a;)  is  an  m  x  n  matrix.  In  this  case,  Qq{x)  =  [Go (a:),  —Im\-  This 
leads  to  the  following  commutative  diagram: 
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The  upper  left  block  is  the  original  problem.  The  lower  right  block  is 
the  new  problem.  This  applies  to  the  eigenproblem,  where  F{x)  =  F(a:„) 
H{x)  -  H(xn)  are  both  square  matrices  that  are  full  rank  except  when  is 
an  eigenvalue. 

Let  Qoix)  =  [A/„, -flirt]  and  show  that  det[F{xn)]  =  0  iflF  det[H{x„)]  =  0. 
This  can  then  be  used  to  prove  that  the  Incomplete  Intersection  algorithm 
ends  up  with  a  new  eigenproblem  whose  roots  are  identical  with  those  of  the 
original  system. 
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